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Preface 


Celestial  mechanics  problems  have  motivated  the  development  of  both  pure  and 
applied  mathematics  since  Newton’s  development  of  calculus.  In  fact,  Newton  used  his 
calculus  to  formulate  Kepler’s  laws  of  planetary  motion  and  to  show  that  a  point  mass 
may  approximate  a  spherically  symmetric  body  when  considering  gravitational  attraction. 
Since  then,  generations  of  mathematicians  have  endeavored  to  improve  the  accuracy  of 
the  equations  for  planetary  motion.  In  doing  so,  analytical  methods,  approximations, 
and  gimmicks  were  formulated,  greatly  influencing  instruction  in  mathematics.  Although 
these  nineteenth  century  mathematicians  could  have  used  the  brute  force  method  of  literal 
numerical  calculation,  without  the  aid  of  modern-day  computers  the  task  would  span 
lifetimes.  Instead,  particular  tricks  found  useful  in  solving  Kepler’s  laws  were  regarded  as 
significant  and  subsequently  added  to  the  mathematician’s  box  of  tools.  In  retrospect  it  can 
be  seen  that  the  precomputer  emphasis  on  analytical  solutions  was  not  entirely  beneficial; 
many  of  these  methods  continue  to  be  taught  as  an  established  mathematical  means  to 
an  end,  despite  an  era  of  sophisticated  computers  where  approximate  methods  have  little 
relevance  (12:421,432),  (6:85).  On  the  other  hand,  analytical  methods  and  their  series- 
2md  closed-form  solutions  yield  insight  into  the  physical  effects  of  specific  parameters  and 
the  relative  strengths  of  various  physical  mechanisms.  Only  the  most  successful  numerical 
computations  provide  physical  insight. 

This  thesis  explores  the  implications  of  a  set  of  mathematical  tricks  employed  for 
centuries  by  those  investigating  celestial  mechanics.  It  appears  that  this  is  one  instance  in 
which  an  analytical  method  produced  a  satisfactory  and  meaningful  result  in  the  past  but 
may  have  also  hidden  the  existence  of  a  great  number  of  alternative  solutions.  However,  the 
spurious  nature  of  these  numerical  results  indicates  that  we  must  not  completely  abandon 
our  centuries-old  mathematical  tricks.  The  emergence  of  supercomputers  has  seemingly 
produced  answers  in  minutes  rather  than  lifetimes,  but  only  slightly  diminishes  our  reliance 
on  mathematical  approximations  and  analytical  methods. 


Cynthia  Ann  Provost 
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Abstract 

Numerous  studies  have  been  conducted  on  equilibrium  orientations  of  objects  moving 
under  the  influence  of  a  central  gravitational  field.  The  objects  under  investigation  have 
had  configurations  ranging  from  spheres  to  cylinders  to  more  complicated  three-axis  stabi¬ 
lized  satellites.  The  results  of  many  of  these  studies  conclude  that  equilibrium  conditions 
exist  only  when  one  of  the  principal  axes  coincides  with  the  radius  vector.  Furthermore, 
these  results  assume  that  the  center  of  force  is  located  within  the  orbit  plane,  thereby  trac¬ 
ing  a  great  circle  orbit.  While  these  previous  works  have  approximated  the  gravitational 
potential,  this  study  examines  relative  equilibria  obtained  by  retaining  an  exact  expression 
for  the  potential  of  a  spherical  primary  body,  as  shown  in  a  recent  paper  by  Wang,  Mad- 
docks,  and  Krishnaprasad.  The  exact  dynamic  equations  for  the  motion  of  a  finite  rigid 
body  in  an  inverse  square  gravitational  force  field  are  investigated.  Only  circular  orbits  for 
a  specific  satellite  model  consisting  of  six  masses  connected  by  three  massless  rigid  rods  are 
considered.  The  system  dynamics  are  comprised  of  seven  nonlinear  equations,  which  were 
numerically  solved  on  a  Cray  computer.  The  existence  of  equilibrium  orientations  which 
establish  non-great  circle  orbits  was  verified  and  other  interesting  results  were  noted.  The 
operational  significance  of  these  results  was  also  examined. 
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STEADY  MOTIONS  OF  RIGID  BODY  SATELLITES 


IN  A  CENTRAL  GRAVITATIONAL  FIELD 

/.  Introduction 

1.1  Background 

The  study  of  motions  due  to  an  inverse  square  gravitational  potential  field  has  been 
explored  for  several  centuries.  Orbiting  bodies  investigated  include  point  masses,  spheres, 
cylinders,  rods,  ball-and-socket  connected  objects,  and  more  complicated  bodies  with  mov¬ 
ing  appendages  and  momentum  wheels.  Until  recently,  many  of  these  studies  have  utilized 
approximate  solutions  to  the  equations  of  motion  obtained  with  a  spherical  primary  body. 
Typically  the  gravitational  force  and  moment  acting  on  the  body  of  interest  has  been 
approximated  with  a  binomial  series  expansion  and  the  center  of  force  is  assumed  to  be 
located  within  the  orbit  plane.  Using  these  approximations  generally  leads  to  a  result 
in  which  equilibrium  conditions  exist  only  when  one  of  the  principal  axes  coincides  with 
the  radius  vector.  This  type  of  orientation  is  known  as  a  gravity-gradient  attitude.  The 
trajectory  for  these  types  of  satellites  is  a  great  circle  orbit. 

A  gravity-gradient-stabilized  satellite  utilizes  the  gravitational  force  of  a  massive 
primary  body  for  stability.  The  gravitational  torques  tend  to  align  the  axis  of  minimum 
moment  of  inertia  with  the  radius  vector  from  the  center  of  mass  of  the  satellite  to  the 
center  of  mass  of  the  primary  body.  Use  of  this  stability  method  requires  that  the  satellite 
have  specific  mass  and  inertia  properties.  The  most  common  feature  of  many  gravity- 
gradient  stabilized  satellites  is  the  presence  of  long  booms.  Some  sketches  of  several  of 
these  types  of  spacecraft  are  depicted  in  Figure  1. 

A  great  circle  orbit  is  one  in  which  the  radius  vector  from  the  center  of  mass  of  the 
primary  body  to  the  center  of  mass  of  the  orbiting  body  traces  out  a  disk.  This  type 
of  orbit  is  depicted  in  Figure  2.  Conversely,  the  radius  vector  of  a  non-great  circle  orbit 
traces  out  a  cone,  as  seen  in  Figure  3.  In  other  words,  “the  center  of  relative  equilibrium 
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ATS-A 


GGTS 


Figure  1.  Sketches  of  Gravity-Gradient  Satellites:  Applications  Technology  Satellite- A 
(ATS-A)  (24),  Gravity-Gradient  Test  Satellite  (GGTS)  (24),  Magnetically  An¬ 
chored  Gravity  System  (MAGS)  (24),  and  the  Naval  Research  Laboratory  - 
164  (NRL-164)(5). 


rotation  does  not  coincide  with  the  center  of  the  force  field  (41).”  For  a  great  circle  orbit, 
it  is  obvious  that  the  radius  vector,  A,  must  be  perpendicular  to  the  angular  velocity 
vector,  fl,  or  X  •  n  =  0.  Further  analysis  reveals  that  both  fl  and  X  must  be  aligned  with 
principal  axes.  However,  this  restriction  is  not  present  in  non-great  circle  orbits.  In  fact, 
it  has  been  shown  that  although  the  offset  from  a  great  circle  orbit  is  very  small  in  certain 
circumstances,  the  attitude  change  associated  with  that  offset  can  be  quite  significant  (42). 

In  a  recent  paper  by  Wang,  Maddocks,  and  Krishnaprasad  (42)  the  concept  of  non¬ 
great  circle  orbits  was  explored.  In  their  study,  the  Hamiltonian  form  of  the  exact  dynamic 
equations  for  the  motion  of  a  finite  rigid  body  in  an  inverse  square  gravitational  force 
field  was  examined.  The  equations  were  found  to  have  a  noncanonical  structure  with 
the  gravitational  potential  as  one  term.  While  previous  works  have  approximated  the 
gravitational  potential,  the  Wang  paper  explores  relative  equilibria  obtained  by  retaining 
an  exact  expression  for  the  potential  of  a  spherical  primary  body.  The  results  of  the 
Wang  investigation  were  promising,  yet  not  definitive.  Using  continuation  techniques, 
nume’  '.s  non-great  circle  orientations  were  discovered  for  an  ideal,  hypothetical  structure 
consisting  of  6  mass  particles.  Seven  nonlinear  algebraic  equations  were  solved  by  successive 
application  of  Newton- Raphson  iterations.  However,  several  discontinuities  in  their  results 
indicated  potential  bifurcation  amd  limit  points,  implying  the  existence  of  distinct  branches 
of  equilibria.  Additionally,  the  operational  significance  and  feasibility  of  non-great  circle 
orbit  applications  was  not  addressed.  Furthermore,  their  results  were  based  on  a  model 
which  may  not  realistically  represent  a  useful  satellite  in  orbit. 

If  these  additional  questions  were  addressed,  the  results  could  have  far-reaching  ef¬ 
fects  on  the  design  of  earth-orbiting  satellites.  Different  dynamical  systems  could  then  be 
analyzed  with  the  exact  potential,  yielding  previously  unknown  equilibrium  orientations. 
Modifications  to  control  systems  which  take  into  account  newly-found  stable  orientations 
might  achieve  better  performance  and  reduce  consumption  of  on-orbit  propellants  used 
for  station- keeping  and  attitude  adjustment.  Perhaps  complex  orbital  missions  could  be 
achieved  with  less  complex  satellite  configurations.  Clearly,  the  results  of  (42)  exhibit  po¬ 
tential  for  future  satellite  design.  However,  further  investigation  is  required  to  more  closely 
examine  the  anomalies  and  significance  of  (42). 
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Figure  2.  A  Great  Circle  Orbit. 


Figure  3.  A  Non-Great  Circle  Orbit. 


1.2  Problem  Definition 

This  investigation  attempted  to  validate  the  findings  of  (42)  and  explore  in  greater 
detail  the  discontinuities  present  in  their  results.  A  computer  program  different  from  that 
used  in  (42),  AUT086  (4),  was  used  to  independently  determine  the  existence  of  non¬ 
great  circle  relative  equilibria  for  the  six-particle  model.  AUT086,  referred  to  as  AUTO 
for  the  remainder  of  this  paper,  is  a  continuation/bifurcation  package,  and  its  use  was 
intended  to  facilitate  the  discovery  of  bifurcation  points  and  successfully  trace  out  branches 
of  relative  equilibria.  AUTO  was  chosen  because  of  recommendations  from  other  users, 
the  existence  of  several  references  and  theses  which  utilized  AUTO,  and  the  software’s 
availability.  An  objective  examination  of  the  operational  utility  of  these  relative  equilibria 
was  also  conducted. 

l.S  Scope 

Although  the  original  intent  was  to  examine  in  further  depth  all  the  results  obtained 
in  (42),  as  well  as  utilize  a  more  realistic  satellite  model,  extensive  difficulties  in  achieving 
valid  and  economically  attainable  results  led  to  a  restriction  in  scope.  Consequently,  this 
study  is  limited  to  the  analysis  of  the  exact  six-particle  molecule  model  used  in  (42). 
As  described  earlier,  only  relative  equilibria  for  the  molecule’s  motion  due  to  an  inverse- 
square  gravitational  field  were  sought.  In  the  continuation  method,  only  one  parameter 
was  varied  in  order  to  keep  the  Jacobian  symmetric,  thereby  eliminating  the  need  for  a 
double  precision  eigenvalue  solver  and  simplifying  computations. 

1.4  Assumptions 

The  dynamical  system  examined  was  assumed  to  have  a  circular  orbit  where  the 
center  of  attraction  may  or  may  not  be  the  center  of  the  orbit.  A  restricted  two-body 
approach  was  used  by  assuming  the  primary  body  is  stationary.  The  solution  obtained 
with  the  earth  as  a  spherically  symmetrical  finite  rigid  body  is  identical  to  the  solution 
obtained  by  treating  it  as  a  point  mass.  Thus,  the  earth  was  modelled  here  as  a  point  mass. 
Consequently,  harmonics  associated  with  the  oblateness  of  the  earth  are  not  considered. 
Although  the  very  small  size  of  the  satellite  model  relative  to  its  orbit  radius  could  lead  to 
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its  treatment  as  a  point  mass,  doing  so  would  prevent  analysis  of  satellite  attitude.  Thus 
the  non-great  circle  relative  equilibria  of  interest  here  were  investigated  by  considering  the 
finite  extent  of  the  orbiting  rigid  body. 


1.5  Method 

The  dynamical  equations  were  numerically  analyzed  on  a  computer,  modelling  the 
rigid  body  satellite  as  a  simple  molecular  structure  with  unequal  moments  of  inertia.  First, 
the  equations  were  coded  and  validated  in  double  precision  on  a  Sun  SparcStation.  Values 
of  the  function  and  the  Jacobian  were  verified  using  the  results  of  the  program  used  in  (42). 
The  code  was  then  incorporated  into  AUTO  subroutines  and  tested  for  proper  integration 
within  the  software.  Subroutines  were  added  to  compute  the  necessary  variables  for  accu¬ 
racy  testing.  Insufficient  accuracy  was  obtained  using  the  Sun  workstation,  therefore  the 
code  was  transferred  to  the  Ohio  Supercomputer  Cray  where  double  precision  accuracy  to 
29  digits  was  possible  (29).  Initial  equilibrium  points  were  obtained  from  values  depicted 
in  Table  2  of  (42).  The  AUTO  code  was  run  numerous  times  with  different  values  for 
stepsizes,  iterations,  and  tolerances  to  achieve  results  with  sufficient  accursu:y.  Hundreds 
of  combinations  were  used  in  an  intensive  search  for  valid  results.  Apparent  numerical 
instabilities  caused  several  problems  in  obtaining  data  and  excessive  computer  processing 
time.  Several  attempts  at  scaling  and  other  modifications  to  AUTO  were  made  to  improve 
computer  efficiency  and  accuracy.  When  this  was  accomplished,  the  results  were  compared 
to  those  found  in  (42).  While  the  results  of  this  project  were  not  as  conclusive  as  had  been 
anticipated,  some  results  of  (42)  were  verified. 

1.6  Overview 

Following  this  introduction  is  a  literature  review  which  gives  a  synopsis  of  signifi¬ 
cant  studies  on  relative  equilibria  of  bodies  moving  under  a  gravitational  field.  Next,  a 
development  of  the  gravity-gradient  orientation  concept  is  presented  in  Chapter  III.  Corre¬ 
spondingly,  a  development  of  the  theory  on  which  this  study  is  based  is  described.  Chapter 
IV  outlines  the  solution  approach  while  in  Chapter  V  the  numerical  results  are  analyzed 
and  interpreted.  Conclusions  r^arding  the  significance  and  operational  applicability  of 
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the  concepts  are  contained  in  Chapter  VI.  Finally,  recommendations  for  furthn  study  are 
provided. 


//.  Littrature  Review 


In  classical  rigid-body  dynamics,  relative  equilibria,  a  term  attributed  to  Pmncard, 
are  stationary  rotations  about  the  principal  axes  of  inertia  (36),  (41).  Although  the  subject 
of  this  report  is  relative  equilibria  of  rigid  body  satellites,  the  concept  parallds  so  closdy 
with  gravity-gradient  satellites  that  a  literature  review  of  the  latter  serves  as  an  excdlent 
background  for  the  former.  Thus,  this  chapter  primarily  focusses  on  research,  operations, 
and  ideas  pertsdning  to  gravity-gradient  torques  and  stabilization. 

2.1  Gravitational  Effects 

According  to  Griffin,  “planetary  gravitational  fields  decrease  with  distance  R  from 
the  center  of  the  planet  according  to  the  Newtonian  1/R^  law,”  assuming  higher  order 
harmonics  are  neglected  (6:293).  Only  the  mass  center  of  the  spacecraft  is  in  a  true 
gravitationally-derived  orbit.  This  explsdns  why  an  orbiting  object  experiences  a  stronger 
attraction  on  its  “lower”  side  than  its  “upper  side.”  A  spacecraft  mass  particle  closer  to  the 
earth  would  move  ahead  of  a  mass  particle  farther  away  from  earth’s  center  if  in  a  free  orbit. 
If  this  differential  force  is  applied  to  a  body  with  unequal  principal  moments  of  inertia,  a 
resulting  torque  tends  to  rotate  the  object,  verticsdly  aligning  its  minimum  inertia  or  longest 
axis.  Alternatively  stated,  “perturbations  from  the  equilibrium  produce  a  restoring  torque 
toward  the  stable  vertical  position,  causing  a  periodic  oscillatory  or  ‘librational’  motion 
. . .  and  internal  elastic  forces  in  the  structure  balance  the  orbital  dynamic  accelerations 
tending  to  separate  masses  orbiting  at  different  altitudes  (6:64-65,293).” 

Several  authors  describe  gravity-gradient  torque  in  a  similar  manner.  Satellites  with 
unequal  principal  moments  of  inertia  have  torques  acting  upon  them  due  to  the  gradient  of 
the  earth’s  gravitational  ii  'if’.  According  to  (39),  “these  torques  tend  to  align  the  axis  of 
least  moment  of  inertia  with  the  local  vertical.  Therefore,  in  principle,  the  gravity-gradient 
torques  can  be  used  as  the  basis  of  a  Siaoilization  system  for  an  earth-pointing  satellite.” 
Thomson  states  that  “gravity  torque  tends  to  align  the  axis  of  minimum  moment  of  inertia 
of  a  non-spinning  body  with  the  radius  vector  drawn  from  the  center  of  E^rth  (38).”  More 
specific  descriptions  of  the  vehicle  orientation  can  be  found  in  (41):  “the  angular  velocity 
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lies  along  the  principal  axis  of  the  body  with  the  largest  associated  moment  of  inertia 
(minor  axis  of  the  ellipsoid  of  inertia),  and  the  radius  vector  is  aligned  to  the  principal  axis 
with  the  least  associated  moment  of  inertia  (major  axis  of  the  ellipsoid  of  inertia),”  and 
(17):  ‘‘maintaining  a  fixed  orientation  with  respect  to  a  rotating  reference  frame  established 
by  the  binormal,  tangent,  and  principal  normal  of  the  orbital  path.” 

The  effects  of  gravity-gradient  torques  can  be  significant,  depending  on  the  distance 
from  the  center  of  the  gravitational  field  and  the  shape  of  the  orbiting  body.  For  large 
vehicles  such  as  a  space  station  or  space  shuttle  the  magnitude  may  be  10~^p’s  or  more 
(6:64).  According  to  (47),  when  considering  gravitational,  magnetic,  solar  radiation,  aero¬ 
dynamic,  and  cosmic  dust  disturbances,  gravity-gradient  torques  are  the  most  dominant 
torque  over  a  large  range  of  altitudes.  Furthermore,  the  spacecraft  yaw  angle  does  not 
influence  the  gravity-gradient  torque,  since  yaw  represents  rotation  around  the  local  verti¬ 
cal.  The  torque  magnitude  depends  on  the  difiference  between  principal  moments  of  inertia. 
Therefore,  cylindrical  spacecraft  are  more  affected  than  disk-shaped  spacecraft  (6:294). 

2.2  Other  Spacecraft  Torques 

To  put  the  gravity-gradient  concept  into  perspective,  a  brief  discussion  of  a  satellite’s 
surrounding  space  environment  is  in  order.  The  review  will  be  limited  to  a  description  of 
various  torques  which  affect  bodies  orbiting  the  earth. 

A  space  vehicle  orbiting  the  earth,  or  any  other  gravitational  body,  is  subject  to 
a  number  of  disturbance  torques.  Although  some  have  postulated  that  the  gravitational 
gradient  torque  is  the  only  external  torque  on  the  body  (17),  there  are  others.  Torques 
due  to  aerodynamic  loads,  magnetic  fields,  solar  wind  and  radiation,  particle  impingement, 
outgassing,  and  internal  moving  parts  also  affect  a  space  vehicle. 

Aerodynamic  torque  arises  from  the  separation  between  the  mass  center  and  aero¬ 
dynamic  center  of  pressure  (6:292).  Momentum  is  transferred  from  atmospheric  molecules 
to  the  surface  of  the  space  vehicle  (13:250).  This  torque  is  also  referred  to  as  atmospheric 
resistance  (30). 
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Magnetic  torque  affects  the  space  vehicle  while  under  the  influence  of  the  earth’s 
magnetic  field.  Comparable  to  the  alignment  of  a  compass  needle  with  the  local  magnetic 
field,  components  aboard  the  spacecraft  which  contribute  to  the  vehicle’s  magnetic  moment 
may  similarly  cause  alignment  with  the  earth’s  magnetic  field.  Consequently,  careful  rela¬ 
tive  orientation  of  current  loops,  magnets,  and  electronics  is  required  by  system  designers 
to  ensure  the  total  m^netic  moment  of  the  vehicle  is  zero  (6:295),  (13:265). 

Solar  rsidiatioa  pressure  torque  also  contributes  to  the  overall  motion  of  a  spacecraft. 
This  torque,  sometimes  attributed  to  the  solar  wind  (31:83),  is  independent  of  velocity  or 
position  while  the  spacecraft  remains  in  sunlight,  and  acts  perpendicular  to  the  sun  line 
(6:295).  At  altitudes  above  1000  km,  the  space  vehicle  disturbance  torque  environment  is 
usually  dominated  by  solar  radiation  pressure  (6:294).  According  to  (49),  solar  radiation 
produces  a  transverse  thermal  gradient  in  deployed  booms,  causing  them  to  bend  away 
from  the  sun.  The  bending  affects  satellites  by  moving  both  the  center  of  mass  and  the 
center  of  solar  pressure,  and  changing  the  inertia  dyadic  which  produces  a  change  in  gravity 
restoring-torque  as  well  as  dynamic  effects  due  to  rigidity. 

Torques  may  also  be  caused  by  passage  through  meteor  showers  or  a  collision  with 
debris.  The  impact  of  a  meteoroid  would  impart  a  short  impulse  to  the  spacecraft  (31:83). 
Additionally,  the  spacecraft  itself  is  the  cause  of  several  disturbance  torques.  Deliberate 
or  accidental  effluent  venting  and  internal  torques  resulting  from  a  momentum  exchange 
also  cause  disturbances.  The  internal  torques  miy  arise  from  the  motion  of  antennas,  solar 
arrays,  instruments,  scanners,  deployable  booms,  and  other  appendages  (6:296).  Those 
appendages  in  turn  would  be  subject  to  linear  acceleration,  centrifugal  force,  and  Coriolis 
force  (30). 

A  graphical  representation  of  some  of  these  torques  can  be  found  in  Figure  4,  which 
originated  in  (47)  and  was  modified  in  (13:271).  The  figure  presents  a  useful  comparison 
of  the  magnitudes  of  some  disturbance  torques  as  a  function  of  altitude.  The  torques 
are  linear  in  the  parameter  shown,  with  the  gravity-gradient  linear  limit  being  about  15°. 
Thus,  the  gravity-  gradient  torque  for  6  =  5*  may  be  obtained  by  taking  the  appropriate 
value  from  Figure  4  and  multiplying  by  5,  where  6  indicates  the  offset  angle  from  the 
vertical. 
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Figoie  4.  Comparison  of  several  common  torques  on  a  typical  spacecraft  (47). 
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2.S  StabUity 

According  to  (24),  the  first  published  concept  outlining  a  passive  gravity-gradient 
control  system  is  believed  to  exist  in  a  U.S.  patent  issued  to  R.E.  Roberson  and  J.V. 
Breakwell  in  1956.  This  patent  was  only  the  first  in  a  long  line  of  technical  reports  con¬ 
cerning  gravity -gradient  stabilization  systems.  One  aspect  of  gravity-gradient  stabilization, 
and  passive  attitude  control  in  general,  is  that  it  ‘‘takes  advantage  of  basic  physical  princi¬ 
ples  and  naturally  occurring  forces  by  designing  the  spacecraft  so  as  to  enhance  the  effect 
of  one  force  while  reducing  others  (6:297).”  Dependable  spacecraft  stabilization  permits 
electronic  and  optical  tracking  system  performance  to  be  optimized  using  more  efficient 
antenna  design  and  optimal  orientation  of  onboard  light  sources  or  refiectors  (18).  In  (1), 
the  author  reduces  the  problem  of  gravity-gradient  stabilization  to  four  major  areas;  1) 
the  reduction  of  tumbling  2)  the  extension  of  a  boom  or  booms  to  achieve  a  large  ratio  of 
pitch-  and  roll-axis  to  yaw-axis  moment  of  inertia  3)  the  damping  of  satellite  librations  to 
tolerable  oscillations,  and  4)  controlling  the  satellite’s  boom  position.  Numerous  authors 
have  SMldressed  these  problems;  several  are  described  in  this  section. 

In  (19),  the  energy  equation  is  used  to  determine  the  stability  extrema  of  asymmetric 
satellites.  Concurring  with  several  other  works,  the  author  states  that  almost  all  satellites 
have  a  stable  orientation.  In  particular,  he  addresses  six  dynamically  distinct  sets  of 
equilibria,  one  of  which  is  stable  for  any  asymmetric  rigid  object  in  a  circular  orbit. 

The  problem  of  relative  equilibria  for  gyrostats  in  a  circular  orbit  is  investigated 
by  (48).  In  his  paper,  the  author  discusses  a  satellite  composed  of  an  arbitrary  number 
of  gyrostats  interconnected  by  ball-and-socket  joints.  A  gyrostat  is  a  stabilizing  device 
consisting  of  a  small  axisymmetric  rotor,  or  momentum  wheel,  spinning  about  its  axis  of 
symmetry  inside  a  platform.  The  rotor’s  spin  does  not  affect  the  mass  distribution  smd 
the  inertia  tensor  of  the  gyrostat  remains  constant  if  the  rotor  and  the  platform  are  rigid 
bodies.  A  gyrostat  is  also  referred  to  as  a  dual-spin  satellite  (31:264).  First  and  second 
variations  of  the  dynamic  potential  energy  are  used  in  (48)  to  obtmn  relative  equilibrium 
positions  and  stability  criteria. 
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The  unpredictability  of  the  direction  of  the  pointing  axis  is  described  in  (43).  The 
author  states:  “All  gravity-gradient  stabilized  satellites  have  a  stable  orientation  which  is 
either  right  side  up  or  upside  down.  If  after  deployment  the  attitude  and  angular  velocity 
are  sufficient  to  cause  tumbling,  the  final  orientation  may  be  upside  down.”  The  energy 
approach  is  used  to  predict  capture  or  tumbling,  thus  indicating  whether  a  turnover  device 
is  required  (43). 

The  lack  of  inherent  yaw  stability  present  with  gravity-gradient  attitude  control  is 
briefly  discussed  in  (6:300).  Such  a  spacecraft  is  unconstrained  and  may  rotate  about  its 
vertical  axis.  Consequently,  large  amplitude  librations  are  sometimes  observed  which  may 
be  sufficient  to  invert  the  spacecraft.  This  anomalous  behavior  has  been  associated  with 
the  long-period  resonances  in  the  gravity-gradient  boom  when  excited  by  thermal  input, 
as  described  in  (5)  and  (28:16). 

In  (38),  the  author  discusses  the  potential  stabilization  of  a  particular  type  of  space¬ 
craft  through  the  use  of  spin.  The  particular  vehicle  configuration,  with  its  symmetric  axis 
normal  to  the  orbit  plane,  is  apparently  unstable  when  /j  <  /].  Ii  and  are  the  moments 
of  inertia  along  the  transverse  and  symmetric  axes  respectively.  According  to  (38),  the 
amount  of  spin  required  to  stabilize  the  spacecraft  depends  on  the  angular  velocity  of  the 
orbit  radius  and  the  ratio  /3//1. 

Presenting  a  general  approach  to  the  “rigorous  nonlinear  stability  analysis  of  rela¬ 
tive  equilibria  in  Hamiltonian  systems”  was  the  author’s  intent  in  (36).  In  this  paper,  the 
author  applied  the  reduced  energy-momentum  method  to  acquire  an  explicit  and  easily  im- 
plementable  criterion  for  nonlinear  stability  of  relative  equilibria.  His  method  involves  only 
the  configuration  variables  and  does  not  introduce  Lagrange  multipliers.  The  method  also 
does  not  require  explicit  knowledge  of  momenta  and  momenta  variations  or  the  conserved 
quantities  in  the  reduced  space,  also  known  as  Casimirs  (36). 

Several  other  authors  have  examined  the  design  of  gravity- gradient  stabilized  space¬ 
craft.  Specifically,  in  (13:281-346),  an  excellent  presentation  on  torques,  stabiUty,  design 
and  flight  experience  pertaining  to  gravitational  forces  is  found.  The  stability  conditions 
obtained  through  variational  methods  is  the  also  the  subject  of  (41).  Additionally,  the 
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effect  of  a  variable  gravitational  force  acting  on  a  rod  satellite  is  presented  in  (31:83). 
Three  types  of  gravity-gradient  stabilization  configurations  —  traac,  magnetic  anchor,  and 
lunged  multibody  —  are  discussed  in  (9). 

While  the  gravitational  torques  present  in  orbit  can  provide  a  useful  stabilizing  force, 
additional  measures  are  sometimes  required  to  put  the  spacecraft  in  a  proper  gravity- 
gradient  orientation  from  its  deployment  attitude,  and  also  to  minimize  the  naturally 
occurring  libration.  Although  it  is  often  required  to  position  a  satellite  in  orbit  oriented  in 
a  specified  direction,  maintaining  the  specified  attitude  may  be  impeded  due  to  perturbing 
torques  (38).  In  (23),  several  methods  for  this  type  of  stabilization  are  described.  Magnetic 
stabilization  can  be  used  to  assist  the  space  vehicle  in  assuming  its  desired  gravity-gradient 
attitude.  The  same  can  be  accomplished  using  the  previously  mentioned  boom  deployment. 
Other  methods  include  magnetic  hysteresis  rods,  a  damping  spring  with  magnetic  hysteresis 
rods,  and  magnetic  hysteresis  plus  eddy-current  damping. 

In  (39),  “a  passive  three-axis  stabilization  system  is  described  which  relies  solely  on 
the  gravity-gradient  to  provide  both  the  damping  and  the  restoring  torques  for  a  vertically- 
oriented  satellite.”  Damping  occurs  through  the  energy  dissipation  of  relative  motion 
existing  between  the  main  spacecraft  vehicle  and  an  auxiliary  body. 

Exploring  allowable  attitude  deviations  due  to  given  disturbances  was  the  intent  of 
the  author  in  (19).  This  information  is  of  particular  use  to  spacecraft  design  engineers  who 
assess  long  term  attitude  performance  of  satellites. 

2.4  Advantages  of  Gravity-Gradient  Stabilization 

There  are  several  advantages  to  utilizing  the  naturally  occurring  gravitational  torque 
in  the  design  of  spacecraft.  A  gravity-gradient  stabilized  spacecraft  will  maintain  a  stable 
orientation  relative  to  the  earth,  without  decaying  or  drifting  due  to  other  environmental 
torques.  This  assumes  a  stable  environment  (45:19).  Gravity-gradient  stabilization  is  one 
of  several  passive  methods  of  stabilization.  Passive  methods  require  no  spacecraft  power 
consumption  or  ground-directed  commands  (45:19).  With  a  position-controlled  satellite, 
mission  requirements  may  be  met  with  a  minimum  of  on-orbit  satellites  and  the  least 
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extensive  and  expensive  system  of  ground  stations  (10).  Another  advantage  of  passive 
attitude  control  is  the  potential  for  very  long  satellite  lifetimes.  The  operational  life  is 
not  limited  by  onboard  consummables  or  wear  and  tear  on  moving  parts  (6:297).  Gravity- 
gradient  stabilization  is  useful  when  long  orbit  lifetimes  are  needed  and  the  attitude  sta¬ 
bility  requirements  are  relatively  broad  (6:300).  The  operational  utility  of  these  satellites 
is  explored  further  in  Section  2.8. 

2.5  Disadvantages  of  Gmvity-Gradient  Stabilization 

While  gravity-gradient  stabilization  does  have  its  advantages,  the  method  is  not 
without  some  detrimental  characteristics.  The  attitude  of  such  a  satellite  is  limited  to  only 
one  or  two  possible  orientations  -  the  minimum  inertia  axis  pointing  toward  or  away  from 
the  center  of  the  gravity  field.  This  method  is  only  effective  near  a  massive  central  body, 
limiting  its  utility  and  preventing  its  application  to  certain  planetary  probes.  Only  specific 
satellite  configurations  permit  the  use  of  this  stabilization  method,  since  it  requires  the 
use  of  long  booms  or  an  elongated  mass  distribution.  A  gravity-gradient  satellite  is  also 
subject  to  wobble  or  libration  (45:19).  Libration  amplitudes  of  10-20®  are  not  uncommon 
(6:300).  Furthermore,  attitude  control  is  limited  to  about  1”.  Thermal  gradients  may  easily 
develop  across  the  length  of  the  satellite  or  its  deployable  booms,  contributing  to  both  the 
libration  and  a  decrease  in  attitude  control  (45:19).  This  problem  is  most  likely  related  to 
the  “unpredictable  low  frequency  instability  associated  with  solar  radiation”  addressed  in 
(5).  In  one  author’s  opinion,  gravity-gradient  stabilization  is  too  imprecise  and  inflexible 
for  most  applications  due  to  its  poor  overall  accuracy  and  somewhat  inflexible  response  to 
changing  environmental  conditions  (6:300). 

2.6  Equilibrium  Orientations 

It  is  at  this  point  that  the  literature  review  becomes  specifically  pertinent  to  this  par¬ 
ticular  study.  Given  the  existence  of  relative  equilibrium  orientations,  one  would  naturally 
speculate  exactly  how  many  equilibrium  attitudes  are  produced  by  gravitational  effects. 
Although  no  articles  professing  the  existence  of  less  than  24  equilibrium  attitudes  were 
located,  several  papers  declared  the  existence  of  24  or  more  equilibrium  attitudes.  The 
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reports  described  in  this  section  all  discuss  whether  or  not  gravitational  effects  produce 
exactly  24  or  at  least  24  equilibrium  attitudes.  This  report  and  (42)  appear  to  resume  the 
debate. 

According  to  (17),  Lagrange  was  apparently  the  first  to  deduce  from  the  relevant 
equations  of  motion  the  existence  of  24  equilibrium  attitudes.  In  these  attitudes  the  body 
principal  axes  are  aligned  with  the  orbiting  reference  frame.  He  also  investigated  satellite 
librational  motions  and  formulated  stability  conditions  corresponding  to  the  linearized 
dynamical  equations.  A  similar  result  is  obtained  in  (19)  where  the  author  claims  that  ‘*a 
rigid  body  without  an  axis  of  symmetry,  in  a  circular  orbit,  possesses  24  equilibria . . .  which 
may  be  classified  into  six  dynamically  distinct  groups  each  containing  four  elements.” 

The  author  of  (17)  continues  by  stating,  “Most  investigators  appear  to  take  for 
granted  the  uniqueness  of  the  24  classical  equilibrium  orientations.  A  contrary  result, 
however,  is  claimed  by  a  more  recent  paper  ...(22)  ...specifically,  the  existence  of  an 
infinite  number  of  attitudes  in  which  an  arbitrary  satellite  can  remain  fixed  in  the  orbiting 
frame.”  The  work  in  (22)  was  refuted  in  (14)  and  again  in  (17)  by  the  same  author.  Another 
paper,  (33),  considers  a  similar  satellite  and  clmms  that  “for  a  solid  body  with  unequal 
moments  of  inertia  moving  in  a  circular  orbit,  there  are  no  equilibrium  positions  other  than 
those  corresponding  to  the  coincidence  of  the  three  principal  central  axes  of  inertia  of  the 
body  with  the  axes  of  the  orbital  coordinate  system.”  His  solutions  correspond  to  different 
cases  for  coincidence  of  the  principal  inertia  axes  of  the  body  with  the  orbital  coordinate 
system  axes.  The  author  then  states  that  “this  result  demonstrates  that  the  findings  of 
Michelson  [(22)]  are  not  correct.” 

In  their  second  attempt  at  describing  the  uniqueness  of  24  equilibrium  orientations, 
the  authors  of  (14)  and  (17)  derive  equations  for  the  attitude  motion  of  a  rigid  body 
satellite  in  a  circular  orbit  about  a  spherically  symmetric  earth.  Additionally,  the  authors 
investigate  the  existence  of  equilibrium  orientations  in  which  the  body  maintains  a  fixed 
orientation  with  respect  to  an  orbiting  frame  of  reference,  orientations  corresponding  to 
gravity-gradient  stabilization  positions  (21).  Intending  to  satisfy  even  the  most  cynical 
non-believers  they  reported  on  a  mathematical  development  of  the  problem  which  reduces 
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to  two  eigenvalue  problems  with  exactly  24  solutions.  Clearly,  not  all  were  convinced,  as 
can  be  seen  in  (40),  (41),  and  (42). 

As  is  described  in  this  paper,  the  authors  of  (40),  (41),  and  (42)  examined  the  equa¬ 
tions  of  motion  pertaining  to  a  rigid  body  of  finite  extent  moving  in  a  central  gravitational 
field.  These  studies  were  intended  to  investigate  the  approximations  inherent  in  models 
of  a  moving  point  mass.  The  authors  questioned  the  common  practice  of  retmning  only 
the  lowest-order,  nonvanishing  term  of  the  gravitational  potential.  While  non-great  circle 
motions  were  shown  to  be  impossible  in  the  point  mass  model,  relative  equilibria  were  ob¬ 
tained  through  careful  consideration  of  the  exact  gravitational  potential  acting  on  a  finite 
body.  The  results  of  these  papers  decidedly  conflict  with  those  of  (14),  (33),  and  (17)  and 
apparently  may  lend  credence  to  the  work  in  (22). 

2.1  UrAque  Ideas 

While  a  coauthor  of  (17)  did  adhere  to  the  notion  of  only  24  equilibrium  attitudes, 
he  also  recognized  the  possibility  of  what  is  referred  to  in  (42)  as  “non-great  circle  orbits.” 
He  admits  that  the  true  orbit  plane  may  be  displaced  from  the  nomiiial  orbit  plane  such 
that  the  true  orbit  plane  does  not  contain  the  gravitating  center  (32).  He  elaborates: 

A  point  mass  in  an  inverse  square  gravitational  field  can  move  in  a  circular  orbit 
whose  plane  contains  the  gravitating  center.  In  fact,  the  gravitating  center  must 
lie  in  the  plane  of  the  orbit,  whatever  conic  the  latter  may  be,  because  there 
is  no  force  mechanism  to  generate  or  sustain  any  out-of-plane  motion.  On  the 
other  hand,  a  real  material  body  is  not  a  point  mass  and  the  gravitational  force 
on  such  a  distributed  body  depends  on  the  orientation  of  the  body  about  its 
center  of  mass.  It  is  quite  possible  for  the  orbiting  body  to  experience  a  force 
which  is  not  coUinear  with  the  line  from  the  attracting  center  to  the  center  of 
mass  of  the  body.  (32) 

As  in  (42),  this  author  indirectly  examined  the  effects  of  approximations  of  the  gravitational 
potential  and  the  different  results  obtained  using  finite  body  rather  than  point  mass  models. 

Finally,  one  other  work  may  have  served  as  a  precursor  to  the  formulation  of  more 
exact  equations  of  motion  for  a  rigid  body  in  a  gravitational  field.  In  (21),  the  effect  of 
higher-order  inertia  integrals  is  investigated: 
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Under  the  assumption  that  the  body  dimensions  are  small  relative  to  the  radial 
distance  from  the  center  of  force,  terms  in  the  third-  and  higher-order  inertial 
integrals  are  ignored.  However  for  a  body  with  all  three  moments  of  inertia 
equal  second-order  terms  in  the  gravitational  torques  are  zero  ...  An  analy¬ 
sis  including  the  higher  order  inertia  integrals  shows  that  the  attitude  of  an 
earth-pointing  satellite  with  three  equal  moments  of  inertia  may  be  rendered 
stable  by  the  differential-gravity  effects,  provided  the  satellite  is  not  spherically 
symmetric  and  no  torques  from  other  sources  exist.  (21) 

Using  the  Liapounov  direct  method  an  expression  for  the  gravitational  potential  was  de¬ 
rived  which  included  inertia  integral  terms  through  the  fourth-order.  His  analysis  showed 
that  approximating  the  gravitational  potential  with  only  second-order  inertia  integrals  is 
not  always  adequate.  This  would  be  particularly  true  when  the  three  moments  of  inertia 
are  equal  or  nearly  equal  (21). 


2.8  Operations 

Utilization  of  gravity-gradient  stabilization  has  been  documented  in  numerous  articles 
and  demonstrated  in  many  fully-functioning  satellites  as  weU.  This  section  describes  a  few 
examples. 

In  June  1963  the  first  orbiting  vehicle  to  successfully  achieve  passive  gravity-gradient 
attitude  stabilization  was  launched.  The  space  vehicle,  designated  1963-22A  or  TRANSIT, 
was  followed  by  1963-38B,  1963-49B,  and  1964-26A  (23).  These  five  spacecraft,  designed 
by  the  Applied  Physics  Laboratory  at  The  Johns  Hopkins  University,  were  instrumental 
in  establishing  the  effectiveness  of  gravity-gradient  stabilization.  Although  one  began 
tumbling  and  stabilized  upside  down,  improvements  were  made  and  lessons  were  learned. 
Stabilization  problems  were  corrected  on  subsequent  missions  with  the  use  of  magnets, 
while  thermal  bending  of  the  100-foot  booms  was  minimized  by  the  use  of  an  alternative 
material  which  was  less  absorptive.  Using  the  silver-plated  booms,  the  12“  amplitude  of 
the  high-frequency,  dynamic  boom  bending  was  decreased  to  less  than  0.25“  (23). 

The  subject  of  booms  arises  frequently  in  the  literature.  In  (5),  it  is  postulated  that 
the  clever  use  of  tip-weighted  extendable  booms  on  a  spacecraft  could  provide  three-axis 
passive  stabilization,  and  thus  a  potentially  desirable  earth-pointing  equilibrium  attitude. 
Gravity-gradient  stabilization  achieved  with  dumbbell  or  boom  configurations  is  also  dis- 
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cussed  in  (31:83).  The  gravity-gradient  loading  on  the  long  booms,  exhibited  during  the 
boom-bending  observed  on  mission  1963-22A,  are  mentioned  also  in  (31:59).  The  torques 
on  the  booms  cause  the  spacecraft  to  assume  a  radial  attitude.  In  (6:299),  the  author 
explains  that  a  typical  method  of  obtaining  the  desired  inertia  properties  for  a  spacecraft 
is  to  deploy  a  ‘^motor-driven  boom  with  a  relatively  heavy  end  mass.”  This  can  be  accom¬ 
plished  by  unrolling  a  reel  of  “prestressed  metallic  tape  —  like  a  carpenter’s  measuring 
tape”  that  forms  a  cylinder  with  significant  lateral  stiffness  but  minimal  torsional  rigid¬ 
ity  (6:299).  This  particular  type  of  boom  is  discussed  in  detail  in  (30).  Since  satellites 
equipped  with  long  rigid  booms  cannot  be  deployed  from  the  fairing  of  a  rocket  nose  cone, 
the  appendages  must  be  extendable  and  retractable  (31:83).  One  solution  to  this  problem 
is  the  STEM  -  a  storable  tubular  extendable  member  (30).  A  STEM  is  a  strip  of  flat,  thin 
material  which  forms  a  very  strong  tubular  shape  when  extended.  Coiled  on  a  drum  when 
stored,  the  STEM  can  provide  a  significant  restoring  torque  as  a  gravity-gradient  stabilizer 
(30). 

Gravity-gradient  stabilized  satellites  are  useful  for  situations  in  which  a  particular 
instrument  is  designed  to  point  in  the  zenith  or  nadir  direction.  This  typically  involves 
a  specific  inertia  configuration  in  which  the  vehicle  possesses  an  axis  where  f,  <<  Is,Iy 
(6:299).  For  example,  the  TRANSIT  radio  navigation  satellite  system’s  main  requirement 
was  a  nadir-pointed  antenna.  These  satellites  exhibited  15-year  lifetimes  (6:297),  (23). 
Although  the  Naval  Research  Laboratory  (NRL)  originally  designed  satellites  to  achieve 
pointing  accuracy  of  ±10“  (1),  more  recent  techniques  have  permitted  accuracies  on  the 
order  of  1-3“  (6:309).  GEOSAT,  a  US  Navy  radar  altimetry  satellite  launched  in  1984, 
achieved  vertical  stability  to  within  1“  using  an  extremely  stiff  boom  with  an  eddy  current 
damper  as  its  top  mass  (6:300).  The  use  of  gravity-gradient  stabilization  on  lenticular 
communications  satellites  is  presented  in  (11). 

Even  if  a  satellite  is  not  designed  to  employ  the  stabilizing  effects  of  gravity-gradient 
torques,  its  effect  on  operations  must  still  be  considered.  For  some  applications,  gravity- 
gradient  effects  may  not  be  of  significance.  However,  certain  industrial  manufacturing 
and  materials-processing  operations  are  best  conducted  in  low  gravity  conditions.  These 
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activities  would  need  to  take  place  near  the  center  of  mass  or  at  a  highm*  altitude  where 
gravitational  effects  are  diminished  (6:64-65). 

To  summarize,  volumes  of  material  have  been  written  on  the  subject  of  gravity- 
gradient  stabilization  and  relative  equilibria  of  rigid  bodies  in  a  central  gravitational  fidd. 
This  chapter  has  provided  a  small  sampling  of  some  of  the  available,  and  possibly  more 
interesting,  results  achieved  in  the  last  three  decades.  Three  particularly  useful  sources  are 
(13),  (28),  and  (24).  For  further  review,  the  reader  is  directed  to  the  bibliographies  found 
in  the  above  references. 
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in.  System  Dynamics 


S.l  Introduction 

First,  a  development  of  the  familiar  gravity-gradient  concept  will  be  ^ven,  makii^ 
special  note  of  the  approximations  made  in  obtaining  a  result.  Next  an  explanation  of  how 
this  investigation  differs  from  the  typical  gravity-gradient  devdopment  will  follow.  Finally, 
the  specific  dynamical  equations  of  interest  to  this  study  will  be  summarized. 


S.2  Gravity-Gradient  Development 

Consider  a  rigid  body,  B,  orbiting  a  large  mass,  B„  as  shown  in  Figure  5.  The 
gravicational  attraction  force  experienced  by  each  particle  of  the  orbiting  rigid  body  is 


GM  R-l-r  . 

-5^ - 5 — —dm 

\R  +  f\^  IR  +  r) 


(1) 


where  (R  +  r)/|R  +  r]  is  the  unit  vector  in  the  direction  from  the  center  of  the  large  body 
to  the  arbitrary  mass  particle,  G  is  the  universal  gravitational  constant,  Af  is  the  mass  of 
the  large  attracting  body,  R  is  the  vector  from  the  center  of  mass  of  the  attracting  body 
to  the  center  of  mass  of  the  orbiting  body,  and  r  is  the  vector  from  the  center  of  mass  of 
the  orbiting  body  to  each  individual  mass  particle  (8).  Simplifying  and  integrating  over 
the  entire  body  to  obtain  the  total  gravitational  force  imparted  onto  the  orbiting  body  we 
see  that 

F  =  -GM  f  ?  ---dm  =  mi.  (2) 

where  indicates  the  integral  over  the  orbiting  ri^d  body. 

As  is  frequently  done  in  astronautical  calculations,  algebraic  manipulation  followed 
by  an  approximation  and  subsequent  Taylor  series  expsmsion  simplify  matters  significantly. 
Algebraic  manipulation  of  the  denominator  of  the  integrand  in  Equation  2  yields 


Figure  5.  A  ri^d  body  orbiting  a  massive  central  body  under  the  influence  of  a  gravita¬ 
tional  field. 


1 


_ 1 _ 

1 


(3) 


Now  we  assume  |r|  «  |R|,  neglect  terms  including  |^^/|R|’,  and  apply  a  Taylor  series 
expansion,  (1  -I-  *)“"  =  1  —  n®,  to  obtain 


|R+fl* 


«  |R| 


-3 


1-34' 


|R|> 


Substituting  this  expression  into  Equation  2  we  get 


/(R  +  r) 
|R|»  Jb^  ' 


1-3?;' 


mi’j 


dm 


(4) 


(5) 
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and 


(6) 


^  GM  .  GM  f 


GM  f  (ft  +  r)(3fi  f) 


|R|» 


lip 


dm 


As  has  been  done  for  centuries,  an  approximation  is  then  made  by  neglecting  higher  order 
terms  and  retaining  only  the  first  nonzero  t«m.  After  acknowledging  that  the  int^al  of 
r  dm  over  the  body  vanishes  we  then  obtain 


F  «  - 


GMm^ 

W 


(7) 


and  finally 

mi  +  ^S„6  (8) 

|R|> 

The  same  result  may  also  be  obtained  if  the  orbiting  body  is  spherically  symmetrical,  or 
originally  defined  as  a  point  mass. 

A  similar  development  of  the  accompanying  moment  equation  and  an  examination 
of  the  balance  of  angular  momentum  for  a  circular  orbit  yield  the  familiar  gravity-gradient 
concept:  by  aligning  two  principal  axes  with  the  orbit  normal  and  the  radius  vector  to  the 
central  body,  equilibrium  is  obtained.  Thus  there  are  24  possible  equilibrium  configurations 
of  a  gravity-gradient  body  in  a  circular  orbit  —  3  axes  to  align  with  the  radios  vector  times 
2  directions  per  axis  times  4  axes  to  align  with  the  orbit  normal  (46),  (17).  Not  all  24 
of  these  are  stable  of  course,  as  is  shown  in  (46:148),  (19),  and  (13:294-307).  Similar 
developments  of  the  gravity-gradient  tensor  can  be  found  in  (45:128),  while  an  expression 
for  total  gravity  gradient  load  as  a  function  of  mass,  length,  angular  velocity  and  vertical 
displacement  angle  is  located  in  (31).  An  expression  for  gravity-gradient  torque  is  also 
given  in  (6:293).  Development  of  the  gravity-gradient  geopotential  and  a  discussion  of 
harmonics  are  provided  in  (20:430-438). 

The  operational  significance  of  this  gravity-gradient  finding  has  been  described  in 
numerous  works  such  as  (28),  (10),  (11),  (23),  and  (6).  Several  satellites  have  incorporated 
this  type  of  passive  stabilization,  namely,  the  Transit-SA,  GGTS  (Gravity  Gradient  Test 
Satellite),  DODGE  (Department  of  Defense  Gravity  Experiment),  GEOS  (Geodetic  Earth 
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Orbiting  Satellite),  RAE^2  (Radio  Astronomy  Ebcplorer),  NRL  164,  Salyut  6,  and  LDEF 
(Long  Duration  Exposure  F^ty)  (45:787),  (5),  (13:327,334,336,339,341). 

S,S  A  Gmvity-Gradient  Alternative 

The  premise  of  this  investigation  is  almost  exclusivdy  based  on  the  results  of  (42). 
This  section  contauns  a  brief  overview  of  the  dynamical  equations  of  interest.  The  reader 
is  referred  to  (40),  (41),  and  (42)  for  a  detauled  account  of  the  theoretical  development  of 
the  relevant  equations  of  equilibrium  discussed  here. 

As  an  alternative  formulation  to  Equation  8,  consider  again  the  motion  of  a  rigid  body 
6*  under  the  influence  of  a  central  inverse-square  gravitational  force  due  to  a  massive  body 
B*.  Because  the  massive  body  is  assumed  to  be  stationary,  a  restricted  two- body  approach 
is  applicable.  An  inertial  reference  frame  is  located  at  the  origin  O,  the  mass  center  of 
the  stationary  massive  body.  As  shown  in  Figure  6,  with  C  designated  as  the  mass  center 
of  the  rigid  body  in  motion,  r*  denotes  the  vector  from  0  to  C  in  the  inertial  frame  and 
Q*  is  the  vector  from  C  to  an  arbitrary  point  in  the  rigid  body  in  the  body  frame.  B  is 
the  transformation  matrix  from  the  body  frame  to  the  inertial  frame.  The  mass  particle 
associated  with  the  arbitrary  point  will  be  denoted  by  dm(Q*). 

It  is  easily  seen  that  the  linear  momentum  and  angular  momentum  of  B  about  O, 
expressed  in  the  inertial  frame  are 

r  =  jf .  +  BQ-)dm(^*)  =  m?*  (9) 

and 

r  =  jr(r*-|-BQ*)x|(r*-|-BQ-)dm(Q*)  (10) 

=  ?•  X  mF*  -I-  Brn*  (11) 

where 
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Figure  6.  A  rifpd  body  in  a  central  gravitational  field. 

m  =  mass  of  the  moving  rigjd  body 

^  s  linear  momentum 

jf*  s  angular  momentum 

n*  =  instantaneous  angular  velocity  in  the  body  frame 

I*  =  the  moment  of  inertia  matrix  of  5*  in  the  body  frame 

and  the  dot  indicates  a  derivative  with  respect  to  time  while  the  superscript  *  denotes  a 
quantity  with  dimension  (for  reasons  which  will  soon  be  apparent).  The  moment  of  inertia 
matrix  I*  is  written  as 

r  =  /  Q-Q-’^dm(Q*)  (12) 

Jb* 

where  Q*  indicates  the  skew-symmetric  matrix 

'  \  /  0  -Ql  Ql  ' 

Qi  =  Ql  0  -Q\  (13) 

\  Qt  0  ^ 
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The  primary  body  is  assumed  to  be  spherically  symmetric  and  the  resulting  potential 
iidd  produced  by  Bl  is 

V=  -  f  — — dm(Q*)  (14) 

Jb-  I?  +  BQ‘|  ^  ^ 

where  M  is  the  mass  of  Bl.  Assuming  that  no  external  forces  are  acting  on  the  two-body 
system,  the  resultant  force  on  5*,  with  P  =  VK,  is 


f  GAf(r* -I- BQ*)  g 


(15) 


This  result  arises  since  balance  of  linear  momentum  implies  that  the  sum  of  forces  is  equal 
to  the  time  rate  of  change  of  linear  momentum.  Note  that  the  final  force  expression  in 
Equation  7  is  different  from  the  one  just  devdioped  in  Equation  15.  Multiplying  this  force 
term  by  a  moment  arm  to  obtain  a  torque,  assuming  there  are  no  external  couples  acting 
on  B*,  yirids 

f=-/  (f*-|-B^*)x  =  d  (16) 

Jb-  ^  ^  |P  -I-  BQ*|»  '  ^  ^ 

which  then  implies,  through  balance  of  angular  momentum. 


if  s:0 


(17) 


Now  the  substitution  n*  =  1*11*  is  made,  with  if*  denoting  the  instantaneous  body 
angular  momentum  of  B*.  Furthermore,  the  equations  are  now  expressed  in  the  body 
frame  through  the  following  vector  body  variables: 


The  equations  of  equilibrium  may  now  be  written  as 


(19) 

‘n*  +  ^ 

m 

(20) 

.  t  GM(V  +  Q*) 

“  L  |V  +  Q-|> 

(21) 

(22) 

The  system  is  then  nondimensionalized  with  mass,  length,  and  time  scales 


=ra‘ 


and  nondimensional  variables 


1  - 

I  ^  “  m/r-‘ 


'4 


By  converting  the  rigid  body  B*  into  a  nondimensional  body  B,  the  following  definitions 
arise  from  the  dimensional  mass  oifferential  dm. 


Q 

di/(Q) 


51 

/ 

ijm(q-) 


lM<i)  =  /  = 

Jb  Jb*  m 


where  the  nondimensional  moment  of  inertia  of  B  is  defined  as 

and  ^r(l)  =  1.  The  dynamic  equations  can  now  be  expressed  with  the  nondimensional 
variables  as 
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(29) 

A  X  r  *n  +  M 

(30) 

(31) 

Bi"»n 

(32) 

with  the  prime  denoting  a  derivative  with  respect  to  the  nondimensional  time  t.  Since 
Equations  29-31  do  not  include  the  attitude  B,  they  are  decoupled  from  Equation  32. 
Also,  since  the  states  of  interest  in  this  investigation  are  those  of  relative  equilibria,  where 
the  rigid  body  B*  spins  steadily  about  O  while  the  center  of  mass  C  remains  in  a  circular 
orbit.  Equations  29-31  may  be  set  equal  to  zero.  Therefore,  the  reduced  equations  of 
equilibrium  are 


n  X  i->n  -h 

/  ^  =  S 

(33) 

Axi-*n-i-M  =  5 

(34) 

I*  X  i-*n  ~ 

(35) 

Through  a  series  of  manipulations  involving  Hamiltonians,  Casimir  functions,  and 
variational  methods,  the  dynamical  equations  are  then  reduced  to 


(I  -I-  A’^A)n 


/3n 


/. 


A  +  Q 

|A  +  Qi> 


d«/(Q) 


c 


(36) 

(37) 

(38) 


where  is  a  li^ange  multiplier  with  a  value  of  —  |A|^  and  c  is  a  constant  from  variational 
principles  (42). 

These  are  the  generalized  equations  of  equilibrium  considered  in  this  report.  Of 
particular  significance  is  the  fact  that  the  gravitational  potential  term  in  these  equations 
was  not  approximated  and  remains  exact.  The  use  of  the  exact  model  is  the  key  difference 
between  this  approach  and  the  gravity-gradient  approach  developed  earlier. 
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As  in  (42),  the  model  used  to  discover  and  demonstrate  these  non-great  circle  orbits 
consists  of  six  mass  particles  connected  by  three  rigid  rods  whose  mass  is  considered  negli¬ 
gible.  Two  of  the  point  masses  lie  at  opposite  ends  of  each  of  the  three  principal  axes.  The 
three  rods  form  an  orthonormal  body  frame  and  coincide  with  the  three  principal  axes. 
The  molecule  model  is  shown  in  Figure  7.  Transforming  the  integral  in  equation  37  to  a 
summation  yields  the  following  model-specific  equations: 


|n|’A-(n.A)n-.4g']ji^-±y 


0 

0 

0 


(39) 

(40) 

(41) 


where  A  is  a  parameter  introduced  to  minimize  numerical  ill-conditioning.  The  value  of 
A  was  chosen  to  be  0.04|A|^.  A  complete  description  of  the  methodology  used  to  arrive 
at  this  value  is  given  in  (42).  The  preceding  equations  serve  as  the  departure  point  from 
which  the  remaining  analysis  proceeds. 
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IV.  Methodology 


4.1  Model  Description 

As  seen  in  Figure  7,  the  model  used  in  this  investigation  consists  of  6  point  masses 
attafhed  to  the  ends  of  3  mutually  perpendicular  massless  rods.  The  exact  configuration 
of  masses  and  rod  lengths  was  chosen  to  maximize  the  effect  of  higher  order  terms,  ac¬ 
complished  by  specifying  moments  of  inertia  close  to  each  other  and  insuring  the  model  is 
highly  asymmetric.  As  described  in  (42),  the  masses  and  moments  of  inertia  were  chosen 
as  follows: 


mi  =  0.330066  m,  =  0.00330033 
m3  =  0.330033  m^  =  0.00330033 
ms  =  0.33  ms  =  0.00330034 

/i  =  0.3332  /j  =  0.3335  I3  =  0.3333  (42) 

These  choices  satisfy  the  constraints  imposed  by  the  nondimensionalization: 

6  6 

23  mi  =  1 

ia]  <=1 

where  the  second  constraint  arises  from  tr(I)  =  1. 

In  order  to  analyze  numerically  the  equations  of  motion  on  a  computer,  equations 
39-41  were  written  in  a  form  which  could  be  readily  converted  to  Fortran  code.  Similarly, 
the  Jacobian,  of  Equations  39-41  was  required  to  perform  the  continuation 

technique  utilized  in  obtaining  relative  equilibria.  The  Fortran  representation  of  both  the 
function  /  and  the  Jacobian  /'  are  provided  in  Appendix  A. 

The  dynamic  equations  were  coded  and  validated  in  double  precision  (REAL*8) 
on  a  Sun  SparcStation.  Validation  included  a  comparison  of  the  values  of  f{X,Sl,0)  and 
/'(A,  SI,  P)  obtained  using  the  code  in  Appendix  A  and  those  obtained  by  using  the  program 
developed  for  (42).  The  function  code  was  then  incorporated  into  an  AUTO  subroutine 
and  tested  for  proper  integration  within  the  software. 


30 


J^.i  Continuation  Method 

Relative  equilibria  were  found  by  solving  the  system  =  0  using  a  Newton- 

Raphson  method.  The  method  will  be  introduced  here  first  in  scalar  form.  As  described 
in  (26),  (29),  (2:153),  and  (7),  the  Newton-Raphson  method  may  be  derived  from  a  Taylor 
series  expansion,  which  can  be  represented  as 

/(a:,-n)  =  /(*.)  +  f\xi){xi+i  -  Xi)  +  ^-^{Xi+i  -  Xif  (44) 

where  ^  is  somewhere  in  the  interval  between  Xi  and  Xi+i.  An  approximation  may  be  made 
by  truncating  Equation  44  after  the  first  derivative  term: 


/(*<+! )  /(a:.)  +  f'{Xi){xi+i  ~  Xi)  (45) 

Since  we  are  trying  to  solve  for  f{x)  =  0,  f(xi+i)  would  be  equal  to  zero  at  the  intersection 
with  the  X  axis.  This  observation  yields 

0  ~  f(xi)  +  /'(x,)(x<+i  -  Xi)  (46) 

which  can  then  be  solved  for  x^+i  as 

The  exact  method  used  in  this  study  follows  a  similar  development,  although  the 
vector  form  is  slightly  different  than  the  scalar  form  described  above.  Beginning  with  an 
initial  known  solution  /(u„,  tjo)  =  0,  where  u  indicates  the  states  and  r;  is  a  parameter,  a 
new  solution  f{u,rj)  is  desired.  First,  an  Euler  step  is  taken  where  ri  =  t)o  +  At;.  We  now 
seek 

f{uo  +  Au,  Tj„  +  At;)  =  0  (48) 

Using  a  multi- variable  form  of  the  Taylor  series  expansion  we  see  that 
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0  =  /(«,  +  Au,tj,  + Ai;) 

=  /(tto, »?«)  +  »7o)Att  +  /,(«,,  i?,)Ai7  +  higher  order  terms  (49) 

where  /„  is  the  matrix  of  partial  derivatives  of  /  with  respect  to  «,  also  known  as  the 
Jacobian,  J,  and  /,  is  the  matrix  of  partial  derivatives  of  /  with  respect  to  ij.  Setting  the 
initial  solution  equal  to  zero  and  solving  for  An  yields 


An  =  %)A»7  (50) 

which  in  turn  gives  a  first  approximation  for  a  new  value  of  u 

J~‘(«»>  Oo)/,(tt*,  i7o)Ai7  (51) 

Now,  a  Newton  step  is  taken  to  improve  upon  the  most  recent  value  of  u.  Keeping 
the  same  value  of  q,  we  now  desire  /(«„  +  An,  q)  =  0.  As  before, 

/(«n  1 »?)  +  •f  («n  > »?)  A«  +  higher  order  terms  =  0  (52) 

and 

A«  =  -  J'‘(«„,  »7)/(«n,  Tf)  (53) 

which  produces  a  better  estimate  of  u  in  the  form  of 

«n+i  =  «„  -  »/)/(«„,  T})  (54) 

Proceeding  in  this  manner  allows  one  to  vary  the  parameter,  q,  and  obtain  successive 
solutions  to  the  original  function  /.  Using  the  notation  of  the  particular  dynamical  system 
examined  in  this  study,  branches  of  equilibrium  points  could  be  determined  by  starting 
with  an  initial  solution  for  /(A,n,/?,c)  and  incrementing  the  parameter,  c.  The  process 
of  estimating  and  improving  successive  values  of  a  state  vector  may  be  referred  to  as 
pre-conditioning.  The  path  following  method  described  here  and  other  methods  of  pre¬ 
conditioning  are  explained  in  further  detail  in  (44). 
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4.S  Utilizing  AUT086 


A  continuation  software  package  called  AUT086  (4),  written  by  Eusebius  Doedel,  was 
used  to  trace  solution  branches  in  this  study.  AUTO  contains  algorithms  capable  of  solving 
a  system  of  algebraic  equations  /(u,  q)  =  0,  where  u  is  a  set  of  states  and  f/  denotes  one  or 
more  variable  parameters.  The  software  is  also  capable  of  detecting  algebraic  bifurcation 
points  and  limit  or  turning  points  along  a  solution  branch.  The  reader  is  referred  to  (35) 
for  an  explanation  of  bifurcation  and  limit  points. 

An  initial  solution  for  a  particular  value  of  the  parameter,  t),  is  used  as  input  for 
AUTO  to  begin  its  continuation  routine.  From  the  starting  point,  the  parameter  is  varied 
in  a  stepwise  fashion  and  a  new  solution  is  traced  using  a  pseudo-arclength  continuation 
technique  (4).  User-supplied  subroutines  are  also  required  to  calculate  the  function  /(u,  q), 
the  Jacobian  /'(u,}}),  and  to  initialize  certain  program  constants. 

From  the  equations 


-  iij-ifuj-i  +  -  As  =  0 


(55) 

(56) 


where  As  is  the  stepsize  along  the  branch,  the  next  solution  is  computed.  and  are 
preassigned  constants  used  to  indicate  scaling  differences  between  u  and  rj.  The  stepsize 
may  be  fixed  or  adaptive.  For  an  adaptive  stepsize  in  which  the  Newton  iteration  routine 
fails  to  converge,  the  stepsize  is  halved.  For  an  adaptive  stepsize  in  which  the  Newton 
iteration  routine  converges  rapidly,  the  stepsize  is  increased.  However,  the  stepsize  is  not 
permitted  to  proceed  beyond  a  user-selected  minimum  and  maximum.  Convergence  criteria 
are  also  selected  by  the  user.  Convergence  is  reached  if  the  Newton  increments  satisfy 


1  +  M  ” 


and 


lIHloo 

1  +  IHI. 


(57) 


where  t,  and  c„  are  selected  by  the  user  (4). 
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Accuracy  Testing 


Some  of  the  AUTO  subroutines  were  modified  and  others  were  added  to  explicitly 
establish  the  accuracy  of  a  computed  equilibrium  solution.  As  described  in  (42),  (25),  and 

(26) ,  the  Kantorovich  theorem  may  be  used  to  measure  the  error  for  a  solution  and  serve  as 
a  criterion  for  convergence.  Consequently,  additional  calculations  were  made  to  determine 
the  necessary  coefficients  for  the  Kantorovich  theorem.  The  following  constraints  were 
used  in  accepting  an  AUTO-derived  equilibrium  point: 

a<T7  <  i  and  7  <  f 

A 

where 

a  =  2n|A| 

n  =:  number  of  equations  =  7 

a  —  ” 

Pmin 

Pmin  =  minimum  eigenvalue  of  /'(X,  n,)9) 

7  =  lIHloo  =  niax(^u<,*=  1.2,...7) 

and  c  is  a  user-specified  error  bound.  If  these  constraints  are  satisfied,  the  order  of  con¬ 
vergence  is  guaranteed  to  be  at  least  quadratic  (25). 

Insufficient  accuracy  was  obtained  using  the  Sun  workstation,  therefore  the  code 
was  transferred  to  the  Ohio  Supercomputer  Cray  where  double  precision  accuracy  to  29 
digits  was  possible  (29).  Even  after  this  was  accomplished,  accuracy  still  did  not  appear 
to  correlate  with  the  results  achieved  in  (42).  Thus,  a  gaussian  elimination  subroutine  in 
AUTO  was  replaced  with  a  singular  value  decomposition  (SVD)  subroutine  obtained  in 

(27) .  The  SVD  routine  was  thought  to  be  a  more  robust  algorithm  well-suited  for  such  an 
ill-conditioned  problem. 

4.5  Procedure 

Significant  importance  was  given  to  reproducing  Table  2  of  (42).  Once  those  results 
could  be  duplicated  and  verified,  continued  analysis  of  the  anomalous  behavior  at  certain 


(58) 
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points  of  the  equilibrium  branch  could  be^n.  In  the  interest  of  clarity,  I^le  2  of  (42) 
is  shown  here  as  Thble  1  where  the  solution  numbers  have  been  added  to  assist  in  iden¬ 
tification.  The  table  is  segmented  according  to  the  three  presumed  equilibrium  branches 
represented. 


Table  1.  Relative  Equilibria  Illustrating  Large  Orientation  Drifts  (42) 


Solution 

|A| 

e(X) 

flKA) 

l«l 

m 

1 

500 

46.8611 

-17.4627 

100 

-54.7456 

-32.6009 

2 

760 

47.8276 

-17.8208 

152 

-54.7761 

-34.1683 

3 

1000 

48.7091 

-18.1277 

200 

-54.8384 

-35.5845 

4 

3000 

55.3232 

-17.2751 

600 

-20.2429 

38.7129 

5 

6000 

65.6627 

-19.5986 

1200 

-15.6572 

22.9702 

6 

8000 

71.0045 

-22.9513 

1600 

-11.4520 

17.2236 

7 

9000 

73.0742 

-25.3278 

1800 

-9.5049 

15.2639 

8 

12000 

78.0382 

22.8848 

2400 

-7.5808 

-10.2578 

9 

15000 

81.2009 

18.9175 

3000 

-6.4284 

-6.8820 

10 

20000 

84.1137 

13.7814 

4000 

-4.9693 

-3.7331 

11 

30000 

86.5362 

8.1721 

6000 

-3.2604 

-1.4159 

12 

40000 

87.5514 

5.6183 

8000 

-2.3792 

-0.7051 

The  numerical  analysis  began  by  converting  the  data  in  Table  1  to  cartesian  coor¬ 
dinates  and  using  a  specific  solution,  or  solution  number,  as  the  initial  starting  point  for 
AUTO.  In  particular,  efforts  were  made  to  trace  equilibrium  branches  which  encompassed 
data  from  Solutions  1  to  4,  Solutions  4  to  1,  Solutions  3  to  8,  Solutions  8  to  3,  Solutions  8 
to  12,  and  Solutions  12  to  8.  Distinction  is  made  between  solutions  with  increasing  versus 
decreasing  A  since  this  method  was  intended  to  locate  the  bifurcation  points  implied  by 
discontinuities  found  between  Solutions  3  and  4  and  Solutions  7  and  8. 

Certain  variables  in  the  AUTO  software  were  extremely  sensitive  and  required  ad¬ 
justment  for  each  initial  solution  in  order  to  produce  an  equilibrium  branch.  Consequently, 
several  hundred  runs  using  a  trial  and  error  process  were  required  to  obtain  results  resem¬ 
bling  valid  data. 

Those  equilibrium  states  that  were  obtsuned  were  converted  to  spherical  coordinates: 
A  =  (rcos^cos^,  rcos^sin^,  rsin^)  for  comparison  to  the  results  found  in  (42),  and  to 
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Figure  8.  Spherical  Coordinate  System. 


aid  in  graphical  representation.  The  spherical  coordinate  system  used  is  provided  in  Figure 

8. 


36 


V.  RestJta 


5.1  Introduction 

Efforts  to  duplicate  Table  1  met  with  partial  success.  In  gmieral,  it  spears  as  though 
the  results  of  (42)  can  be  validated,  but  not  with  the  ease  expected  from  using  the  AUTO 
software  package.  The  intent  was  to  perform  a  complete  turning  pmnt  and  bifurcation 
analysis  using  AUTO,  particularly  in  the  repons  where  Table  1  indicated  discontinuities. 
A  few  bifurcation  p<unts  were  detected  but  certain  convergence  tests  were  not  met  and  the 
data  was  deemed  invalid.  The  data  obtained  thus  far  indicates  that  the  scdutions  in  'Dible  1 
are  in  fact  equilibrium  solutions,  but  the  continuous  branches  suggested  by  (42)  could  only 
be  reproduced  by  using  the  same  program  employed  in  (42).  Preliminary  results  implied 
the  existence  of  not  three,  but  at  least  twelve  branches  corresponding  to  each  solution  of 
Table  1;  however,  this  outcome  was  attributed  to  numerical  instabilities  present  in  the 
utilization  of  AUTO  for  this  particular  dynamical  system.  All  results,  including  those 
considered  defective,  were  obtained  at  significant  computational  expense,  possibly  due  to 
these  instabilities.  Specific  results  are  described  in  the  following  section. 

5.2  Discussion  of  Results 

As  stated  earlier,  a  primary  objective  of  this  study  was  the  reproduction  of  selected 
results  from  (42),  shown  in  Table  1.  Figure  9  is  a  graphical  representation  of  Table  1 
with  missing  data  supplied  by  interpolation  and  represented  by  dashed  and  dotted  ( — ) 
lines.  Figure  10  depicts  the  results  obtained  from  several  versions  of  code  and  numerous 
program  runs.  Throughout  the  figures  in  this  section,  data  from  Table  1  is  also  included 
for  comparison  purposes.  Furthermore,  of  the  six  major  equilibrium  solution  components, 
A,  9(A),  ^A),  (2,  9(n),  and  ^12),  only  two,  ^fl)  and  A,  were  chosen  to  graphically 
represent  these  results  since  a  plot  of  ^(2)  most  clearly  indicates  the  existence  of  three 
distinct  equilibrium  branches. 

As  can  be  inferred  from  Figure  10,  only  a  partial  mapping  of  the  complete  phase  space 
was  accomplished.  However,  AUTO  was  able  to  converge  on  the  initial  starting  sdutions 
obtained  from  Table  1,  but  quickly  diverged  from  the  expected  solutions.  Many  equilibrium 
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♦(Q). 


Figure  9.  Interpolated  Data  from  Table  1. 
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X- dimensionless 


x10^ 

Figure  10.  The  results  from  AUTO  are  depicted  as  a  thick  line  of  connected  dots  while 
the  interpolated  results  from  Table  1  are  displayed  by  asterisks  and  dash-dot 
lines.  Additional  AUTO  results  are  found  close  to  the  Table  1  cluster  located 
between  -30  and  -40  degrees,  but  are  not  visible  using  this  scale. 


solutions  were  found  which  cannot  be  seen  in  Figure  10  due  to  the  chosen  scale.  These  are 
explored  in  later  figures.  Additionally,  an  equilibrium  branch  not  obtained  in  (42)  appears 
to  exist  between  values  of  -20*  and  -30*.  The  additional  branch  was  obtained  when 
NEWT  failed  to  converge  on  a  nearby  solution.  NEWT  is  the  program  originally  used 
in  (42)  to  iterate  on  an  initial  state  using  Equations  39-41  until  convergence  on  a  nearby 
solution  occurs.  An  AUTO  run  starting  from  Solution  4  terminated  prematurely.  The 
final  solution  obtained  from  this  run  was  used  as  an  initial  state  in  NEWT,  which  then 
converged  on  a  very  different  solution.  It  is  this  state  which  was  used  as  a  starting  solution 
for  the  results  found  between  -20*  and  -30*.  Section  5.5  elaborates. 

The  following  presentation  of  results  expands  upon  Figure  10  by  dividing  the  region 
into  five  ranges  corresponding  to  areas  of  interest  in  Table  1.  Range  1  includes  equilibria 
found  with  A  values  between  500  and  1000,  corresponding  to  Solutions  1  and  3.  Similarly, 
the  remaining  ranges  are  allotted  according  to  Table  2. 


Table  2.  Definition  of  ranges  used  in  displaying  results. 


Range 

Min  A  value 

Max  A  value 

Starting  Solution 

Ending  Solution 

1 

500 

1000 

1 

3 

2 

1000 

3 

4 

3 

9000 

4 

7 

4 

12000 

7 

8 

5 

12000 

40000 

8 

12 

Figure  11  shows  the  equilibria  obtained  with  Solutions  1,  2,  and  3.  Figure  12  shows 
an  enlargement  of  the  equilibria  obtained  with  an  initial  starting  solution  from  Solution 
1.  The  solutions  near  -32.65*  were  obtained  after  AUTO’s  restart  option  did  not  function. 
The  last  available  solution  was  then  iterated  upon  in  NEWT  and  the  new  starting  solution 
it  produced  was  used  for  an  additional  AUTO  run.  It  should  be  noted  that  although  the 
AUTO  results  appear  to  produce  a  constant  the  variable  does  decrease  as  expected, 
although  not  at  the  rate  implied  by  interpolated  Table  1  results.  AUTO  output  with 
Solution  3  as  the  starting  solution  is  shown  in  Figure  13. 

Results  from  the  end  of  Range  2  are  shown  in  Figure  14.  The  starting  solution  for 
this  data  was  Solution  4.  In  an  effort  to  explore  the  discontinuity  present  in  the  data, 
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Figure  11. 


^12)  -  degrees 


X- dimensionless 


Figure  12.  Enlarged  AUTO  results  in  Ran^  1  are  depicted  by  dots  while  interpolated 
data  from  Table  1  is  shown  by  a  dash-dot  line. 
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Figure  13.  AUTO  results  in  Range  2  are  depicted  by  dots  while  interpolated  data  from 
Table  1  is  shown  by  a  dash-dot  line. 


each  solution  obtained  from  this  run  was  iterated  upon  in  NEWT,  thereby  arriving  at  a 
NEWT-reiined  solution  from  each  AUTO  solution.  The  results  of  the  refinement  process 
are  shown  in  Figure  15.  The  NEWT  data  has  been  partially  filtered  to  eliminate  spurious 
results  attributed  to  the  lack  of  significant  digits  available  on  a  Sun  workstation.  As  would 
be  expected,  the  solutions  are  nearly  identical.  However,  it  is  interesting  that  the  slope  of 
the  NEWT  data  is  of  higher  magnitude  than  the  AUTO  data. 

Figure  16  depicts  the  AUTO  solutions  shown  in  Figure  14  and  those  obtained  from 
a  NEWT  starting  solution.  As  mentioned  above,  after  an  unsuccessful  restart  attempt  the 
last  available  AUTO  solution,  with  Solution  4  as  the  initial  state,  was  used  in  NEWT  to 
arrive  at  a  very  different  solution.  Executing  AUTO  with  this  starting  solution  yielded  the 
results  shown  in  Figure  16.  Figure  17  is  an  enlarged  view  of  the  potential  new  branch(es). 
Although  Figure  17  depicts  four  potential  branches,  the  discontinuities  are  assumed  to  be 
a  result  of  numerical  instabilities  identical  to  those  which  produced  the  results  in  Figure 

14.  Furthermore,  the  difference  between  the  lower  two  branches  and  upper  two  branches  of 
Figure  17  can  be  attributed  to  the  restart  and  NEWT-refinement  process  described  above 
and  displayed  in  Figure  12.  Consequently,  the  four  branches  in  Figure  17  may  represent 
the  potential  existence  of  only  one  equilibrium  branch. 

The  results  obtained  for  Range  4  are  depicted  in  Figure  18.  Again,  it  is  clear  that 
using  AUTO  for  this  particular  dynamical  system  does  not  yield  expected  results.  The 
numerical  instabilities  inherent  to  this  problem  appear  to  cause  AUTO  to  diverge  from 
the  true  equilibrium  solutions.  Again,  it  should  be  noted  that  the  AUTO  results  axe 
not  constant,  they  only  fail  to  vary  at  the  same  rate  as  in  (42).  The  results  shown  in 
Figure  19  were  obtained  in  a  similaur  manner  to  those  of  Figure  12,  with  the  discontinuities 
produced  by  the  restart  and  NEWT-refinement  process  described  earlier.  Eaxh  solution 
composing  the  line  at  -10.2°  was  used  ais  input  for  NEWT  amd  a  refined  set  of  data  was 
obtained.  While  there  were  mamy  spurious  solutions,  just  as  described  with  regard  to  Figure 

15,  the  overwhelming  majority  of  NEWT-relined  solutions  compared  favorably  with  the 
interpolated  Table  1  data. 

Accepting  the  values  obtained  in  (42)  as  valid  equilibrium  solutions,  one  can  see  that 
a  number  of  equilibrium  attitude  orientations  may  now  be  possible.  Figure  20  depicts 
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<t>(Q)  -  degrees 


X- dimensionless 

Figure  14.  Additional  AUTO  results  in  Ranges  2  and  3  are  depicted  by  dots  while  inter¬ 
polated  data  from  Table  1  is  shown  by  a  dash-dot  line. 
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-  degrees 


^(Q)  -  degrees 


Figure  17.  Enlarged  view  of  lower  portion  of  Figure  16.  AUTO  results  in  Range  3  are 
depicted  by  dots  while  interpolated  data  from  Table  1  is  shown  by  a  dash-dot 
line. 
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Figure  18.  AUTO  results  in  Range  4  are  depicted  by  dots  while  interpolated  data  from 
Table  1  is  shown  by  a  dash-dot  line. 
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^(£2) -degrees 


Figure  19.  AUTO  results  in  Range  5  are  depicted  by  dots  while  interpolated  data  from 
Table  1  is  shown  by  a  dash-dot  line.  These  solutions  were  obtained  by  starting 
with  point  a,  which  corresponds  to  Solution  8.  Discontinuities  present  at 
points  b  and  c  were  the  result  of  restart  dilRculties  and  the  NEWT-refinement 
process. 
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a  typical  gravity-gradient  orientation,  as  well  as  the  three  genial  attitudes  indicated  in 
Table  1.  The  molecule  model  is  not  to  scale,  but  representative  of  the  asymmetrical  modd 
used  for  this  study.  With  the  center  of  the  gravitational  body  located  vertically  down, 
the  three  non-great  circle  orientations  are  distinctly  different  from  the  typical  gravity- 
gradient  configuration.  These  equilibrium  attitudes  clearly  indicate  new  possibilities  for 
future  satellite  design. 

5.S  Modifications  to  Improve  AUTO 

Modifications  were  made  to  the  AUTO  software  itself  to  enhance  its  performance 
for  this  particular  dynamical  system.  Since  the  Cray  did  not  support  a  complex  double 
precision  eigenvalue  solver,  AUTO  subroutines  were  modified  to  take  advantage  of  the 
symmetric  Jacobian  in  this  system.  The  symmetric  Jacobian  ensured  real  rather  than 
complex  eigenvalues,  thus  a  complex  double  precision  eigenvalue  solver  was  unnecessary. 

Numerical  instabilities  observed  later  in  the  study  prompted  the  use  of  a  singular 
value  decomposition  (SVD)  rather  than  gaussian  elimination  subroutine  to  solve  the  system 
of  equations.  It  was  believed  that  the  robust  SVD  method  would  eliminate  some  of  the 
problems  addressed  in  Section  5.5.  Although  accuracy  was  improved,  as  indicated  by 
the  Kantorovich  variables  in  Equations  58,  the  SVD  subroutine  produced  nearly  identical 
solutions  to  those  obtained  with  the  gaussian  elimination  subroutine. 

In  order  to  achieve  results  as  close  as  possible  to  those  found  in  (42),  an  additional 
convergence  test  was  placed  in  AUTO  to  ensure  the  desired  accuracy  of  solutions.  As 
described  in  Section  4.4,  in  order  to  satisfy  the  Kantorovich  Theorem  conditions,  7,  the 
maximum  value  of  6u,  was  required  to  be  sufficiently  small.  While  AUTO’s  convergence 
tests,  shown  in  Equation  57,  were  being  satisfied  and  implied  accurate  results,  frequently 
the  Kantorovich  conditions  were  not  met.  One  measure  of  effectiveness  observed  was  the 
norm  of  /,  which  should  be  a  very  small  number  since  the  goal  of  the  program  was  to  find 
successive  solutions  to  /(A,  12,/?)  =  0.  The  value  of  fnerm  obtained  with  AUTO  was  not 
as  small  as  that  achieved  by  using  NEWT.  Therefore,  an  additional  restriction  was  added 
to  the  acceptance  criteria  for  AUTO-derived  equilibrium  solutions.  Additional  iterations 
were  the  result  if  7  was  not  sufficiently  small.  As  expected,  processing  time  increased 
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and  smaller  values  (rf  fn»rm  were  obtained,  although  still  not  as  small  as  those  achieved  in 
NEWT. 

5.4  Additional  Steps  Toward  Efficiency  and  Accuracy 

Given  that  supercomputer  time  is  not  without  cost,  several  attempts  were  made 
at  finding  methods  to  reduce  Cray  processing  time.  In  order  to  ascertain  if  use  of  a 
supercomputer  was  absolutely  essential,  a  comparison  was  made  between  the  Cray  results 
and  the  Sun  SparcStation  results.  Unfortunately  the  comparison  proved  the  Cray  would 
indeed  be  necessary  to  obtain  any  worthwhile  results.  Using  the  same  initial  solutions, 
the  values  of  /n«rm  for  one  iteration  on  the  Cray  and  Sun  respectively  were  0.00001  and 
50,000.  Likewise,  the  values  of  007  were  0.01  and  6,424  —  indicating  that  few,  if  amy, 
accurate  solutions  would  be  found  using  the  Son.  The  AUTO  subroutines  which  define  the 
system  of  equations  to  be  scdved  was  made  as  efficient  as  possible,  minimizing  calculations 
and  reducing  compiling  time.  The  subroutine  which  computes  the  value  of  /(A,  fl,  p)  may 
be  called  over  a  thousand  times  per  Euler  step,  thus  economy  of  calculation  is  imperative. 
In  an  effort  to  increase  the  range  of  solutions  obtained  in  a  given  number  of  steps,  the 
variable  A  was  considered  to  be  the  parameter  rather  than  e.  However,  this  approach 
was  abandoned  since  the  associated  Jacobian  would  then  become  asymmetrical,  causing 
significant  problems  with  the  agenvalue  solver  as  described  in  Section  5.3. 

Furthermore,  several  scaling  schemes  were  tried  after  expected  results  were  not  ob¬ 
tained.  It  was  realized  late  in  the  study  that  AUTO  actually  solves  the  8  x  8  matrix  formed 
from  Equations  55  aund  56.  The  prior  assumption  was  that  only  the  7x7  matrix  indicated 
in  Equations  39-41  amd  55  was  analyzed  in  AUTO.  Additionally,  because  the  value  of  u(7) 
or  P  is  not  pairt  of  AUTO’s  output,  it  was  not  initially  obvious  that  the  values  of  tt(l)  - 
u(7)  ranged  significantly  in  magnitude.  For  these  reasons,  efforts  were  made  to  scale  the  7 
initial  values,  then  the  Jacobian,  and  finally  the  value  of  the  function  /(A,  fl,  P).  Attempts 
to  obtain  convergence  with  AUTO  following  the  scaling  modifications  were  unsuccessful. 
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5.5  Problems  Encountered 


The  most  identifiable  problem  was  an  inability  to  obtain  sidutions  over  a  large  range 
in  a  reasonable  amount  of  CPU  time.  In  order  to  acquire  solutions  with  sufficient  accuracy 
in  AUTO,  the  stepsize  had  to  be  made  so  small  that  the  computational  time  required  to  find 
solutions  was  entirely  too  excessive.  For  example,  achieving  a  continuation  from  A  =  100 
to  103.9  required  over  15  minutes  of  CPU  time  on  the  Cray.  Clearly,  pursuing  this  inves¬ 
tigative  route  would  have  been  impractical.  In  order  to  fully  utilize  AUTO’s  capabilities, 
further  ‘iine-tuning”  of  the  equilibrium  equations  in  the  subroutines  is  required. 

Acquiring  a  suitable  stepsize  for  use  in  AUTO  also  presented  significant  problems. 
A  number  of  variables  in  AUTO  may  be  adjusted  by  the  user  to  achieve  varying  degrees  of 
efficiency  and  accuracy:  mesh  sizes;  initial,  minimum,  and  maximum  stepsizes;  maximum 
number  of  newton  iterations;  and  tolerances  —  almost  a  dozen  altogether.  The  results 
were  extremely  sensitive  to  changes  in  these  factors.  Several  hundred  combinations  of 
these  variables  were  used  in  a  trial  and  error  procedure  throughout  the  course  of  this 
study  in  an  attempt  to  find  an  efficient  method  of  achieving  results.  Many  tradeoffs  were 
made  between  the  variables,  as  decreasing  a  tolerance  might  increase  the  accuracy  but 
would  also  prevent  initial  convergence  on  the  second  solution  or  cause  a  CPU  time  limit 
to  be  exceeded. 

Another  idiosyncrasy  observed  during  this  particular  application  of  AUTO  was  dis¬ 
covered  while  attempting  to  increase  the  stepsize.  An  effective  measure  of  success  was 
obtained  by  comparing  the  parameter  increase  in  two  runs  of  the  same  number  of  steps. 
During  this  process  it  was  noted  that  different  results  were  achieved  depending  on  how 
often  the  data  was  printed  to  the  screen.  The  value  of  the  tenth  calculated  solution  ob¬ 
tained  by  executing  20  steps  and  printing  each  one  was  greater  than  the  value  of  the  tenth 
calculated  solution  obtained  by  executing  20  steps  and  printing  every  other  solution.  The 
source  of  this  discrepancy  was  never  determined. 

Several  other  puzzling  issues  arose  while  employing  AUTO.  Although  a  positive  step- 
size  was  supposed  to  cause  an  increase  in  the  parameter  varied  by  AUTO,  the  reverse  was 
true  for  this  system  of  equations.  Mathematically  this  appears  impossible,  but  the  intrica- 
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des  of  this  unique  system  of  equations  combined  with  AUTO’s  algorithms  appear  to  have 
caused  this  anomaly.  Additionally,  although  AUTO  has  the  capability  to  restart  using 
an  initial  value  obtained  from  any  portion  of  a  previously  obtained  equilibrium  branch,  in 
this  case  the  feature  was  successful  only  sporadically.  This  presented  a  significant  challenge 
since  the  numerical  instabilities  in  the  problem  caused  program  execution  to  terminate  pre¬ 
maturely  numerous  times  due  to  computer  processing  time-limit  violations.  The  numerical 
instabilities  themselves  may  also  be  the  cause  of  the  restart  problem.  Without  a  functioning 
restart  capability,  an  alternative  avenue  for  continuing  the  branch  was  required. 

Possible  causes  for  the  restart  problem  were  analyzed.  It  was  speculated  that  the 
initial  solution  provided  to  AUTO  was  not  of  suiiident  accuracy  to  create  a  valid  equilib¬ 
rium  branch  with  which  the  restart  option  would  work.  Consequently,  the  conversion  of 
initial  states  from  Table  1  was  recalculated  in  double  precision  on  the  Cray  and  also  refined 
using  NEWT,  once  it  was  apparent  that  using  double  predsion  did  not  solve  the  problem. 
This  technique  was  moderately  successful  in  achieving  better  results.  However,  the  restart 
problem  persisted  and  the  refinement  process  was  extended  to  include  the  last  equilibrium 
solution  determined  by  AUTO,  as  well  as  the  initial  solution.  The  presumption  was  that 
AUTO  had  strayed  too  far  from  the  actual  equilibrium  branch  and  was  unable  to  converge 
on  a  subsequent  solution.  Occasionally  the  refinement  process  worked,  but  often  NEWT 
was  unable  to  converge  on  a  nearby  solution  —  the  program  either  failed  completely  or 
provided  a  solution  very  different  from  the  initial  equilibrium  solution.  Steps  were  then 
taken  to  retain  a  certain  number  of  AUTO  solutions  in  double  precision  (29  significant  dig¬ 
its)  as  possible  restart  values.  AUTO,  as  written,  only  retains  solutions  with  10  significant 
digits  for  use  in  subsequent  restarts.  Unfortunately,  this  method  was  also  only  occasionally 
effective.  Another  possible  cause  for  the  restart  problem  is  that  the  states  were  extremely 
close  to  one  another  and  the  minuscule  variation  in  those  states  indicated  a  bifurcation 
point  to  AUTO,  preventing  a  successful  restart.  As  mentioned  earlier,  a  likely  explanation 
for  the  restart  problem  concerns  the  numerical  instabilities  associated  with  these  equilib¬ 
rium  equations  and  the  large  number  of  significant  digits  required.  The  most  apparent 
solution  to  this  problem,  keeping  these  results  independent  of  those  of  NEWT  and  (42), 
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would  involve  program  execution  for  hundreds  of  thousands  of  steps  and  require  the  use 
of  batch  files.  Due  to  limited  computer  time,  this  approach  was  regarded  as  unwise. 
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VI.  Conclusions 


6.1  Summary 

Selected  results  from  (42)  were  verified,  thereby  confirming  the  existence  of  non-great 
circle  relative  equilibria  for  a  rigid  body  traveling  in  a  central  gravitational  field.  Solutions 
were  obtained  utilizing  AUT086,  a  continuation/bifurcation  software  package.  Anticipated 
bifurcation  points  and  limit  points  were  not  found,  presumably  due  to  numerical  sensitiv¬ 
ity  and  instabilities  present  when  employing  AUTO  to  analyze  this  particular  dynamical 
system.  Using  AUTO  and  NEWT,  a  similar  program  implemented  in  (42),  values  and 
trends  of  equilibrium  solutions  were  nearly  identical  to  those  obtained  in  (42).  Ebctensive 
execution  of  the  AUTO  program  and  the  current  subroutines  asso  'iated  with  this  system 
could  possibly  duplicate  and  elaborate  upon  the  complete  results  of  (42),  however  this 
was  not  attempted  due  to  prohibitive  computational  processing  time.  It  is  believed  that 
modification  of  the  current  subroutines  may  lead  to  a  more  efficient  approach  and  thus 
more  definitive  results. 

6.2  Operational  Signifiamce 

The  utilization  of  relative  equilibria  in  non-great  circle  orbits  holds  significant  po¬ 
tential  for  both  current  and  future  satellites.  Current  spacecraft  are  subject  to  attitude 
control  restrictions  based  on  the  classic  gravity-gradient  model.  Greater  natural  attitude 
stability  for  those  satellites  with  active  stabilization  systems  occupying  non-great  circle  ori¬ 
entations  can  increase  pointing  stability  and  accuracy.  These  satellites  would  require  less 
artificial  or  active  attitude  control  and  reduce  attitude  perturbations  inherent  in  artificial 
means  (15).  Missions  might  be  expanded  without  suffering  a  significant  loss  in  pointing 
accuracy.  Satellites  employing  a  mass  expulsion  stabilization  system  could  achieve  and  sus¬ 
tain  additional  and  more  complex  attitudes  naturally,  thus  reducing  propellant  required 
for  attitude  control  and  increasing  the  spacecraft’s  operational  lifetime.  Small  space  test 
spacecraft  with  limited  lifetimes  of  days  or  months  might  profit  considerably,  since  these 
types  of  satellites  would  not  likely  have  complex  active  stabilization  systems  onboard. 
They  might  rely  instead  on  an  expanded  range  of  stable  attitudes  caused  by  gravitational 
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torques.  A  three-axis-stabilized  spacecraft  may  have  been  designed  for  one  set  of  equilib¬ 
rium  points  but  may  now  be  modified  to  operate  in  a  new  range  of  equilibrium  points. 
More  importantly,  prior  knowledge  and  incorporation  of  the  additional  range  of  stable  at¬ 
titudes  would  aid  engineers  in  designing  satellites  with  enhanced  performance.  Reduced 
station- keeping  propellant,  greater  maneuverability,  and  longer  lifetimes  could  be  realized. 
Throughout  the  design  phase,  new  attitudes  and  configurations  could  be  considered  in 
planning  experiments  without  sacrificing  payload  requirements. 

The  missions  which  might  employ  non-great  circle  relative  equilibria  range  from 
astronomy,  sp'xe  test  experiments,  and  space-based  radar,  io  remote  sensing,  surveil¬ 
lance,  and  reconnaissance.  Non-great  circle  relative  equilibria  could  be  useful  for  stellar 
mapping  in  the  infrared  and  visible  wavelengths.  The  instrument  could  be  placed  in  a 
sun-synchronous  orbit  with  its  major  axis  pointed  away  from  the  sun  and  the  natural 
gravitational  torque  could  be  used  for  passive  attitude  stabilization  (16).  For  example, 
although  the  11,400-kg  Hubble  Space  Telescope  is  in  a  28.45'*-inclination,  610-km  orbit 
as  a  result  of  shuttle  payload  and  launch  constraints,  a  more  ideal  orbit  would  have  been 
sun-synchronous  (37),  (34),  (3),  (16).  The  current  orbit  allows  stellar  observations  for 
only  a  portion  of  the  orbit  —  that  part  in  which  the  earth  sufficiently  obscures  the  sun’s 
rays  —  considerably  decreasing  useful  observation  time.  A  sun-synchronous  orbit  utilizing 
these  new  relative  equilibria  for  attitude  stabilization  might  permit  continuous  coverage  of 
deep  space.  Of  course,  in  order  for  a  similar  deep  space  telescope  to  reach  this  orbit,  an 
alternative  to  the  shuttle  would  be  required  for  launch  and  deployment. 

Further  examination  of  non-great  circle  relative  equilibria  could  possibly  explain 
anomalous  behavior  of  current  gravity-gradient  satellites,  providing  that  eccentricity,  geopo¬ 
tential,  thermal,  and  atmospheric  effects  have  already  been  eliminated  as  potential  causes. 
A  spacecraft  exhibiting  unexplained  motion  away  from  the  classical  gravity-gradient  atti¬ 
tude  may  actually  be  drifting  toward  these  new  equilibrium  points.  While  it  is  typically 
thought  that  classic  gravity-gradient  orientation  subjects  a  spacecraft  to  a  zero  overall 
gravitational  torque  load,  perhaps  residual  torque  exists  and  the  vehicle  is  inclined  to  shift 
toward  an  attitude  of  zero  load.  Proof  of  this  concept  would  then  necessitate  modiftcation 
to  Figure  4.  While  the  figure  indicates  that  a  satellite  at  25,228  km  (10'^  km)  experiences 
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a  gravitational  torque  of  8.1721  x  10'“  N-m  if  at  an  orientation  8.1721”  from  the  vertical, 
an  orientation  obtained  from  Solution  11  of  Table  1,  the  satellite  may  actually  be  in  a  state 
of  relative  equilibrium  and  thus  experiencing  zero  net  gravitational  torque. 

A  major  benefit  to  non-great  circle  relative  equilibria  is  yaw  control.  As  discussed  in 
Section  2.3.  gravitational  torques  have  little  influence  on  motion  about  the  vertical  or  yaw 
axis.  However,  in  a  non-great  circle  orbit  the  off-axis  directions  of  the  radius  and  angular 
velocity  vectors  prescribe  a  specific  orientation  for  an  asymmetric  satellite,  thus  possibly 
improving  yaw  control. 

Unfortunately,  despite  the  utility  of  non-great  circle  orbits,  it  is  doubtful  that  any 
onboard  satellite  computer  or  ground-generated  command  signals  would  be  capable  of 
maintaining  one  of  the  relative  equilibrium  attitudes  described  in  this  study.  A  current 
estimate  of  attitude  pointing  accuracy  is  placed  at  one  arcsecond  (16).  Translated  into 
significant  digits  of  attitude  coordinates  this  would  equate  to  six  orders  of  magnitude  — 
considerably  fewer  than  the  29  significant  digits  seemingly  required  by  a  supercomputer  to 
locate  the  equilibrium  points.  Clearly,  a  more  accurate  control  system  would  be  necessary 
to  take  advantage  of  these  newly-found  equilibria.  However,  there  is  some  consolation  in 
noting  that  at  present,  this  represents  a  technological  challenge  rather  than  an  obstacle 
imposed  by  the  laws  of  nature. 

6.3  Recommendations 

Undoubtedly,  this  subject  requires  further  investigation.  It  appears  that,  in  this 
particular  case,  employing  AUTO  to  trace  equilibrium  branches  and  locate  bifurcation 
and  limit  points  is  not  an  uncomplicated  task.  Two  potential  courses  of  action  are  rec¬ 
ommended  with  regard  to  AUTO.  First,  a  thorough  investigation  of  the  AUTO  code  for 
identifiable  sensitivities  in  double  precision  is  suggested.  Particular  attention  should  be 
given  to  insuring  that  all  constants  are  defined  in  double  precision.  Closer  examination  of 
the  code  in  (42)  may  be  required  in  order  to  decipher  the  calculations  used  to  determine 
the  value  of  the  function  /(u).  Subtle  differences  in  mathematical  expressions  for  sum  may 
have  lead  to  discrepancies  between  the  codes  used  in  this  study  and  in  (42).  Additional 
examination  of  the  use  of  the  SVD  subroutine  is  also  suggested.  Apparently,  in  NEWT’s 
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application  of  SVD,  the  small  elements  of  the  diagonal  matrix  W,  the  “small  w(j)’s,’’  are 
not  zeroed  as  is  advised  in  (27:55-56).  Further  experimentation  with  changes  in  mesh  pa¬ 
rameters  (NTST  and  NCOL)  is  also  advised.  Most  importantly,  a  suitable  scaling  scheme 
must  be  developed  in  order  to  improve  pr<^ram  efficiency.  Given  finite  computational 
and  economical  resources,  the  scaling  issue  is  critical.  The  satellite  model  could  also  be 
changed  to  reduce  the  asymmetrical  characteristics  which  magnify  the  numerical  instabili¬ 
ties  present  in  the  problem.  Second,  if  the  first  course  of  action  is  futile,  use  of  an  alternate 
continuation  code  is  recommended. 

Once  the  results  of  (42)  are  extensively  v<didated,  whether  through  AUTO  or  an 
original  program,  the  equations  of  equilibrium  should  be  analyzed  for  a  solid  body  which 
more  closely  resembles  an  orbiting  vehicle.  One  potential  candidate  is  the  NRL-164  (5). 
Accomplishment  of  this  objective  would  further  establish  potential  applications  to  real 
satellites  and  be  a  prelude  to  full-scale  development  and  operational  use  of  spacecraft 
employing  non-great  circle  relative  equilibria. 
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Appendix  A.  AUTO  Subroutines 

This  appendix  contains  the  mathematical  expressions  used  in  creating  the  AUTO 
subroutines  for  this  thesis.  The  associated  code  is  also  included. 

The  vector  form  of  Equations  39-41  was  converted  into  scalar  form  as  follows: 


fi  =  (^?  "i"  ^2  "i"  —  (OiAj  +  OjA]  4"  f23A3)f2i  —  A  8Ufn\  (®®) 

f2  —  "i"  flj  +  fl3)A3  —  (fliAj  fl2A2  "i"  fl3A3)03  —  A  sum]  (®0) 

~  (flj  "I"  Oj  -|-  03}A3  —  (f2iA|  -j-  QjA]  f23A3}03  —  A  sung  (^1) 

/»  =  -fiifli  +  (Aj  +  Aj  +  A3)Oi  —  (OiAi  +  djA]  -|-  d3A3)Ai  +  /3ili  (62) 

/s  =  /jid]  +  (Aj  +  Aj  -f-  ADd]  —  (diA]  +  djA]  +  d3A3)A2  fiilg  (63) 

fe  =  Iggilg  +  (Aj  -|-  Aj  -|-  A3)d3  —  (diAi  +  djA]  +  d3A3)A3  +  /iflg  (64) 

fr  =  i(d?  +  d|-|-d|)-c  (65) 


where 


sumi 


+ 

+ 

+ 

+ 

+ 


sum]  = 


+ 

+ 

+ 


_ ”*i[Ai  +  Qi(l)] _ 

[(Ai  +  gi(l))>  +  (A]  +  Qi(2))*  +  (A3  +  <?i(3))»]i 

_ _ ^2[Ai  +  ^2(1)] _ 

l(Ai  +  <32(1))*  +  (A]  + 1?2(2))*  +  (A3  +  Q2(3))*]^ 

_ rnsjAi  +  ^3(1)] _ 

[(Ai  +  ^3(1))^  +  (A]  +  <?3(2))*  +  (A3  +  Q3(3))*]^ 

_ ^^[Ai  +  Q4(1)] _ 

((Ai  +  <94(1))^  +  (A]  +  ^4(2))*  +  (A3  +  Q4(3))*]^ 

_ ”*5[Ai  +  Qs(l)] _ 

[(Ai  +  ^5(1))*  +  (A2  +  ^5(2))*  +  (A3  +  Q6(3))*]^ 

_ ^6[Ai  -f  Q6(1)] _ 

[(Ai  +  Q6(1))^  +  (A,  +  ^6(2))^  +  (A3  +  Q6(3))2]i 

_ ^i[A2  -f  Qi(2)] _ 

[(Ai  +  Qi(l))^  +  (A]  +  Qi{2)y  +  (A3  +  Qi(3))^]^ 

_ ^2[A2  +  ^2(2)] _ 

[(Al  +  Q2(l)r  +  (^^2  +  <?2(2))2  +  (A3  +  <?2(3))2j« 

_ ”*3[A2  +  ^3(2)] _ 

[(Al  +  ^3(1))^  +  (A]  +  ^3(2))*  -I-  (A3  +  ^3(3))*]^ 

_ ”*<[A2  +  ^4(2)] _ 

[(Al  +  Q,{i)y  +  (A]  +  Q,{2)y  +  (A3  +  g4(3))2]i 


(66) 
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»n«[Aa  +  ^5(2)] 

[(Ai  +  <?5(1))’  +  (A,  +  Q,(2)r  +  (Aa  +  <?5(3))»]f 

m«[Aa  +  ^6(2)] 

[(Ai  +  g,(l))>  +  (Aa  +  ge(2))*  +  (A,  +  Qeim^ 

_  ”*|[A3  +  ^1(3)] _ 

[(Ai  +  Q,{l)y  +  (Aa  +  Q,(2))»  +  (A,  +  Qi(3))*]4 

_ *»a[A3  +  <?a(3)] _ 

((Ai  +  ga(l))*  +  (Aa  +  <?a(2))»  +  (A3  +  ga(3))»]i 

,  »Ti3[Aa  +  ^3(3)] _ 

[(Ai  +  ^3(1))’  +  (Aa  +  ^3(2))*  +  (A3  +  g3(3))*]« 

_ ”*4[A3  -f  ^4(3)] _ 

((A,  +  <?,(!))>  +  (Aa  +  g,(2))2  +  (A3  +  ^4(3))*]^ 
_ »»8[A3  +  ^5(3)] _ 

l(Ai  +  ^5(1))’  +  (Aa  +  ^5(2))*  +  (A3  +  g6(3))*]l 

,  0*6{A3  -t-  ^5(3)] _ 

[(Ai  +  ^6(1))’  +  (Aa  +  ^6(2))*  +  (A3  +  ^6(3))’]^ 


(67) 


(68) 


The  Jacobian,  f\X,{l,0),  of  Equations  39-41  was  also  required  to  perforin  the  nu¬ 
merical  analysis  of  the  equations  of  equilibrium.  The  Jacobian  may  be  expressed  as 

{(|n|’-E?=ii:fe)i-nn’'  -(A  n)i  +  2AnT-nA’'  0 

+ EL.  ife(A +«.)(*+<?,)’■} 

-(A  •  0)1  -I-  2An'*'  -  OAT  I  +  (|A|2  -  p)l  -  -n 

0  -no 

(69) 

where  the  vector,  rather  than  component,  form  is  used  in  the  interest  of  simpUcity. 

The  following  pages  contain  a  listing  of  applicable  subroutines  used  in  AUTO  for  this 
particular  dynamical  system.  However,  the  SVD  subroutines  are  not  included  since  they 
are  copyrighted  in  (27). 


/'(A,n,/?) 
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SUBROUTIIE  FUlGClDIN.U.ICP.PiR.IJAC.F.OFDU.DFDP) 


C  . 

c 

C  This  8ubroutin«  avaluatas  the  right  hand  side  of  the  first  order 
C  systMi  and  the  derivatives  with  respect  to  (U(1),U(2)} 

C  and  with  respect  to  the  free  paraaeters. 

C 

C  Input  paraneters  : 

C  HDIM  -  Dimension  of  U  and  F. 

C  U  -  Vector  containing  U. 

C  PAR  -  Array  of  parameters  in  the  differential  equations. 

C  ICP  -  PARdCPd))  is  the  initial  'free*  parameter. 

C  PAR(ICP(2))  is  a  secondary  'free*  parameter, 

C  for  subsequent  2-parameter  continuations. 

C  IJAC  -  ■!  if  the  Jacobians  DFDU  and  DFDP  are  to  be  returned 

C  >0  if  only  F(U,PAR)  is  to  be  returned  in  this  call. 

C 

C  Values  to  be  returned  : 

C  F  -  F(U,PAR)  the  right  hand  side  of  the  ODE. 

C  DFDU  -  The  derivative  (Jacobian)  sith  respect  to  U. 

C  DFDU(i,j)  must  be  given  the  value  of  d  F(i)  /  d  U(j) 

C 

C  DFDP  -  The  derivative  vith  respect  to  the  ’free’  parameters: 

C  DFDP(i,ICP(j))  =  d  F(i)  /d  PAR(ICP(j)). 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CSGLE  IMPLICIT  REAL  (A-H,0-Z) 

C 

COMMON  /ERRAN/  GANNA,  ALANDA,  SIGNAE,  ALPHA,  RH0NIR,FR0RN 
DOUBLE  PRECISION  U(NDIN) ,PAR(20) 

DOUBLE  PRECISION  F(NDIN) ,DFDU(NDIN,NDIM) ,DFDP(NDIN,20) 

DOUBLE  PRECISION  In(3,3) ,m(6) ,mag21Q(6) ,laiB2,num 
DOUBLE  PRECISION  8um(3) ,dsdl(3,3) ,0(6,3) 

INTEGER  I,J 

C  DATA  STATEMENTS 

DATA  ((d8dl(i,j),i»l,3),j*l,3)/9  ♦  0.0/ 

DATA  ((In(i,j) ,j-l,3) ,i»l,3)/. 3332, 0,0,0, .3335,0,0,0, .3333/ 

DATA  ((Q(i,j),j=l,3),i»l,6)/18*0.0/ 

DATA  (m(i) ,i*l,6)/0.330066d+0,0.00330033d+0,0.330033d+0, 

+0 . 00330033d+0 , 0 . 33d+0 , 0 . 00330034d+0/ 
c 

0(1,1)  *  7.0731919616600413122865203958d-2 
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C' 


Q(2,l)  •  -7.0738992101313602961134447941d-K> 

Q(3.2)  •  7.067S314289714873863597724646d-2 
Q(4.2)  -  -7.0675314289714838201090703958d40 
Q(5.3)  -  7.072489203617885186701998299d-2 
Q(6.3)  -  -7.071760597982097S6822S6668867d-K) 

I****  Ua«ful  quantitlas  a**************** 

«ii2  ■  u(4)**2+u(5)**2+u(6)**2 
omdlan  ■  u(4)*u(l)+u(5)*u<2)+u(6)*u(3) 
laffi2«u(l}**2  ♦  u(2)**2  ♦  u(3)**2 
A*0 . 04* (sqrt (laB2) ) **5 

C**********  Calculating  the  sub  ********* 

do  10  j«l,3 
8ua(j)*0.0 
do  20  ial,6 

num«B( i) * (u( j ) +9 ( i , J ) ) 

mag21Q(i)-(u(l)+Q(i,l))**2+(u(2)+Q(i,2))**2  ♦ 
(u(3)+Q(i,3))**2 
dan>8qrt (Bag21Q ( i) **3) 

8Ub( j )«8ub( j )*nuB/den 
20  continue 
10  continue 

C***********  Calculating  F  *************************** 
do  30  j«l,3 

F( j )»om2*u( j ) -OBdlaB*u( j*3) -A*8um( j ) 

F(j+3)»(In(j ,l)*u(4))+<In(j ,2)*u(5))+(In(j ,3)*u(6))+ 
(lam2*u(j+3))  -(oBdlaB*u(j))+(u(7)*u(j+3)) 

30  continue 

F (7) ■ . 5*OB2-par ( 1) 

f norm«F (1) **2+F (2) **2+F (3) **2+F (4) **2+F (5) **2+F(6) **2 
.  ♦F(7)**2 
f  noxB«sqrt  (f  norm) 

C  srite(6,*)  'fnorm*’ ,fnorm 

IF(IJAC.Eq.0)REnJRN 

C 

C  *****  Calculating  the  partiads  of  suaO  with  respect  to 
C  laBbda(l,2,3) 

C 

DO  26  1-1,3 

d8dl(i,i)«d8dl(i,i)  +  (m(j)*(mag21Q(j))**(-l.S))  -  3.0  * 
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■(j)  •  (Mg21Q<j))**(-2.B)*(u(i)>Hi(j,i))**2 

25  Continu* 

d»dl(l,2)-d«<il(1.2)  -3.0*(*(j)*(Mg2m(j))**(-2-6))  * 
(u(l)+Q(j  ,l))*(u(2)*Q(j  ,2)) 
dsdl(l,3)-dadl(1.3)  -3.0*(B(j)«(ug21Q(j))**(-2.5))  * 
(u(l)+Q(j .l))*(tt(3)*Q(j ,3)) 
dsdl(2.3)-dsdl(2.3)  -3.0*(B(j)*(Bag21Q(j))**(-2.5))  * 
(u(2)+Q(j ,2))*(u(3)*Q(j ,3)) 

26  Continua 

d8dl(2.1)-dBdl(l,2) 

dBdl(3.1)-dsdl(1.3) 

dsdl(3,2)-<lsdl(2.3) 


C  D«rivativ«s  of  F()  with  rospoct  to  U() 

C 

DFDU(l.l)-  u(5)**24u(6)**2  -  A*d8dl(l,l) 
DF0U(1.2)-  -u(4)*u(S)  -  A«<Udl(1.2) 
0FDU(1.3)-  -u(4)*u(6)  <  A*dBdl(1.3) 
DFDU(1,4)-  -u(6)*u(2)  -  u(6)*u(3) 
DFDUd.S)-  2.0*u(S)*u(l)  -  u(4)*u(2) 
DFDU(1,6)-  2.0*u(6)*u(l)  -  u(4)*u(3) 
DFDU(1,7)«  0.0 

DFDU(2,1)-  -u(4)*u(5)-A*<Uldl<2.1) 
DFDU(2,2)-  u(4)**2  ♦  u(6)**2  -  A*d«dl(2,2) 
DFDU(2,3)-  -u(5)*u(6)  -  A*dsdl(2,3) 
DFDU<2,4)-  2.0*u(4)*u(2)  -  u(6)*u(l) 
DFDU(2,6)-  -u(4)*u(l)  -  u(6)*u(3) 
DFDU(2,6)-  2.0*u(6)*u(2)  -  u(6)*u(3) 
DFDU(2,7)-  0.0 

DFDU(3,1)-  -u(4)*u(6)  -  A*dsdl(3.1) 
DFDU(3,2)-  -u(S)*u(6)  -  A*d8dl(3.2} 
DFDU(3,3)=  u(4)>»*2  +  u(5)**2  -  A*d8dl(3,3) 
DFDU(3,4)»  2*u(4)*u(3)  -  u(6)*u(l) 
DFDU(3,S)-  2*u(5)*u(3)  -  u(6)*u(2) 
DFDU(3,6)«  -u(4)*u(l)  -  \i(5)*u(2) 
DFDU(3,7)«  0.0 
DFDU(4,1)«  DFDU(1,4) 

DFDU(4,2)«  DFDU(2,4) 

DFDU(4,3)-  DFDU(3,4) 

DFDU(4,4)»  In(l,l)+u(2)**2  ♦  u(3)**2  ♦  u(7) 
DFDU(4,6)«  In(l,2)-u(2)»u(l) 

DFDU(4,6)»  In(l,3)-u(3)*u(l) 

DFDU(4,7)«  u(4) 

DFDU(6,1)-  DFDUCl.S) 

DFDU(5,2)-  DFDU(2,5) 
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DFDU(5.3)« 

DFDU(3.6) 

0FDU(5.4)- 

In(2,l)-u(i)*u(2) 

DFDU(5.5)« 

In(2,2)-*'u(l)**2-Hi(3)**2  ♦  u(7) 

DFDU(S.6)- 

In(2.3)-u(3)*u(2) 

DPDU(5.7)- 

u(5) 

DFD0(6.1)- 

DFDU(1.6) 

DFDU(6.2)« 

0FDU(2.6) 

DFDU(6,3)- 

DF0U(3.6) 

DFDU(6,4)- 

In(3,l)-u(l)*u(3) 

DFDU(6.5)- 

In(3,2)-u(2)*tt(3) 

DFDU(6.6)- 

In(3,3)-^u(l)**2-»u(2)«e2  ♦  u(7) 

DFDU(6,7)- 

u(6) 

DPD0(7,1)- 

0.0 

DFDU<7,2)« 

0.0 

DFDU(7,3)- 

0.0 

DFDU(7,4)- 

u(4) 

DFDU(7,6)» 

u(S) 

DFDU(7,6)- 

u(6) 

DFDU(7,7)» 

0.0 

c 

C  D«rivatiT«s  vith  raspact  to  PAR(l). 
C 

DFDP(1,1)»0.0 
DFDP(2.1)»0.0 
DFDP(3.1)-0.0 
DFDP(4,1)»0.0 
DFDP(5.1)-0.0 
DFDP(6,l)-0.0 
DFDP(7.1)— 1.0 


RETURM 

E>D 


C 


SUBROUTIIE  STPRTCRDIM.U.PAR) 


C  - 

c 

C  In  this  subrout ina  tha  staady  state  starting  point  anist  ba  defined. 

C  (Used  when  not  restarting  from  a  previously  computed  solution) . 

C  The  problem  parameters  (PAR)  may  be  initialized  here  or  else  in  IMIT. 


C 

C 

C 

C 

C 

C 


HDIM  -  Dimension  of  the  system  of  equations. 

U  -  Vector  of  dimension  RDIN. 

Upon  return  U  should  contain  a  steady  state  solution 
corresponding  to  the  values  assigned  to  PAR. 

PAR  Array  of  parameters  in  the  differential  equatioxis. 
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c 

C  u.dat  is  s  6-lins  input  dsts  fils  consisting  of  ths  3  conpon«its 

C  of  lanbds  snd  ths  3  coi^nsnts  of  oasgs  (in  csrtssisn  coordinatss) 

C 

IMPLICIT  DOUBLE  PRECISIQI  (A-H.O-Z) 

CSGLE  IMPLICIT  REAL  (A>H.Q-Z) 

C 

C 

DIMEISIQI  U(HDIM) .PAR(20) 

OPEB  (UMIT-4,  FILE  •  'u.dat*) 

DO  22  1-1.6 
READ  (4,33)  U(I) 
vrits(6,33)  u(i) 

22  COITIHUE 
33  fomat(D36.30) 

U(7)— (u(l)s*2  ♦  u(2)**2  ♦  u(3)*s2) 
par(l)-0.Sd0*(u(4)**2+u(6)**2+u(6)**2) 

CLOSE  (UMIT  «  4) 

44  RETURM 
ERD 
C 

SUBROUTIHE  IMIT 

C  - 

c 

C  In  this  subroutins  ths  ussr  should  sst  thoss  constants  that  rsquirs 
C  valuss  that  diffsr  from  ths  dsfault  saluss  assignsd  in  DFIRIT. 

C  (Sss  ths  auiin  documsntation  for  ths  dsfault  assignasnts) . 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

CSGLE  IMPLICIT  REAL  (A-H.O-Z) 

C 

COMMON  /BLBCN/  NDIN.IPS.IRS.ILP.ICP(20) ,PAR(20) 

COMMON  /BLCDE/  NTST,NCOL,IAD,ISP,ISN,IPLT,NBC,NINT 
COMMON  /BLTHT/  THETAL(20).  THETAU 
COMMON  /BLDLS/  DS.DSMIN.DSMAX.IADS 
COMMON  /BLEPS/  EPSL(20) .EPSU.EPSS 
COMMON  /BLLIM/  NMX.NUZR.RLO.RLl.AO.Al 
COMMON  /BLMAX/  NPR.MXBF.IID.ITNX.ITNH.NUTN, JAC 
C 

HDIM-7 

IPS-1 

ILP-1 

ISP-2 

C 

C  Set  ths  principal  bifurcation  paraastsr  to  bs  PAR(l). 
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c 

ICP(1)«1 

c 

C  S«t  th«  ■•cond  fr—  paraMtar  to  ba  PAR(2). 
C  (For  2-paraMtar  continuation). 

C 

C  ICP(2)»2 

NZBF-SO 

ITia-200 

IlD-0 

ISV-1 

■UZR-0 

OPEN  (UIIT-S.  FILE  -  'in.datO 

READ  (5.*)  IRS 

READ  (S,*)  ITST 

READ  (5.*)  ICOL 

READ  (S.*)  DS 

READ  (5.*)  DSMIl 

READ  (S.*)  DSNAX 

READ  (5,*)  INZ 

READ  (5.*)  ITRH 

READ  (5.*)  IPR 

READ  (5.*)  THETAU 

READ  (S.*)  THETALCl) 

READ  (S.*)  EPSU 
READ  (5,*)  EPSL(l) 

READ  (5,*)  EPSS 
READ  (5.*)  RLO 
READ  (5,*)  RLl 
READ  (5.*)  AO 
READ  (5,*)  A1 
READ  (S,*)  HXBF 
READ  (5,*)  ITHX 
READ  (5.*)  IID 

CLOSE  (UIIT»5) 

RETURN 

END 

C 

FUNCTION  USZR(I.NUZR.PAR) 

C  - 

c 
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O  Ci 


C  This  subroutiss  can  ba  usad  to  obtain  plotting  and  raatart  data 
at  cartain  valuaa  of  fraa  paranatars. 

IHPLICIT  DOUBLE  PRECISIQI  (A-H.O-Z) 

CSGLE  IHPLICIT  REAL  (A-H.O-Z) 

C 

DIMEMSIOH  PAR(20) 

C 

C  Initially,  for  tha  ataady  atata  analyaia.  aat  IUZR»0  in  HIT. 

C  Than  tha  functions  apacifiad  balov  sill  ba  ignorad. 

C 

C  During  tha  sacond  run  of  this  tast  problaa,  whan  coaqputing  tha  branch 
C  of  par iodic  solutions,  sat  RUZRa4  in  HIT. 

C  In  this  axanpla,  output  sill  than  ba  srittan  in  unit  8  for  tha  saluas 
C  of  PAR(ll)  apacifiad  balos. 

C  Rota  that  PAR(li)  is  normally  rasarvad.  It  is  usad  by  AUTO  to  kaap 
C  track  of  tha  pariod  (Saa  main  documantation) . 

C 

0010(1,2,3,4)1 

C 

1  USZR»PAR(11}  -  11.0 

C  RETURN 

C 

2  USZR«PAR(11)  >14.0 

RETURN 

C 

3  USZR»PAR(11)  -  20.0 

RETURN 

n 

t 

4  USZR«PAR(11)  -  30.0 

RETURN 

C 

END 

C  - 

SUBROUTINE  CAB1(U,DU) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z} 

CSGLE  IHPLICIT  REAL  (A>H,0>Z) 

C 

C  daterminas  tha  max  valua  of  DU  (or  dx)  for  Kantorovich  tast 
C 

DIMENSION  DU(7),  U(7) 

COMMON  /ERRAN/  GANNA,  ALANDA,  SIGNAE,  ALPHA,  RH0NIN,FN0RM 
GAMMA«0.0d0 
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DO  10  1-1,7 

IF(DABS(DU(I)).GT.GAIIKA)  GAIOU-DABS(OU(I)) 
10  COHTIIUE 

ALA)IDA»SQRT(U(1)**2  *  0(2} **2  ♦  0(3) **2) 
c  MUTE (6.*)’ GANNA  -  ’.GAMMA 

RETORN 
END 


C  - - - 

SOBROOTIIE  CAB2(EV) 

C 

IMPLICIT  OOOBLE  PRECISIOI  (A-H.O-Z) 

CSGLE  IMPLICIT  REAL  (A-H.O-Z) 

C 

C  dotaninaa  rhoain,  sigaa,  alpha  for  Kantororich  taat 
C 

DIMEMSIOR  EV(7) 

COMMON  /ERRAN/  GAMMA,  ALANDA,  SIGMAE.  ALPHA.  RHONIM.FHORN 

RHOMIH-1.0<120 

DO  10  1-1,7 

IF(ABS(REAL(EV(I))) .LT.RHOMII)  RHONIH-ABS(REAL(EV(I))) 

10  COHTIMOE 

SIGMAE  -  7.0/RH0MII 
ALPHA  -  2.0*  7.0*  ALAMDA 
RETORN 
END 


C  - - - 

SOBROOTINE  CAB3(R0RLM.DRLN ,RL .RDOMX.DONX ,ONX ,EPSL ,EPSO) 
C 


IMPLICIT  DOOBLE  PRECISION  (A-H,0-Z) 

CSGLE  IMPLICIT  REAL  (A-H,0>Z) 

C 

C  compares  values  vith  tolerances  to  determine  convergence 
C 


DIMENSION  RL(20),  EPSL(20) 

IF(RDRLN.LE.EPSL(1) . AND. RDONX.LE.EPSO) then 
write(6,*)  ’(DRLM-  ’ .DABS(DRLN).’)’ 

write(6,*)  » -  -  '.RDRLM,' 

Brite(6,*)  '1+(RL(1)-',DABS(RL(1)),’)’ 
write(6,*)  ’(DOMX  -  *,  DOMX,’)* 


write(6,*)  ' -  -  ',RDOMX,» 

write(6,*)  '1+(0MX-»,0MX,’)’ 


endif 

RETORN 


<-  ’.EPSLd) 

<-’,EPSO 
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o 


EID 


C  - 

SUBROUTIIE  CAB4 
C 

IMPLICIT  DOUBLE  PRECISIOI  (A-H.O>Z) 

CSGLE  IMPLICIT  REAL  (A>H.O-Z) 

C 

C  prints  kantorovich  raluss  just  prior  to  convsrgsncs 
C 

COMMON  /ERRAI/  GAMMA.  ALAMDA.  SIGMAE.  ALPHA,  RHOMIH.FMORM 

«rite(6,*)  *Fnon  ■  fnonn 

HRITE(6.*)  >GAMMA  -  GAMMA 

WRITE  (6.*)  >ALPHA*SIGMA*GAMMA  -  *.  ALPHA*SIGMAE*GAMMA 

IF(ALPHA*SIGMAE*GAMMA.GE.0.5)  THEN 

WRITE(6,«)  'ALPHA*SIGNA*GAMMA  IS  OUT  OF  BOUNDS  ************> 

ENDIF 

RETURN 

END 


SUBROUTINE  CAB5(NN,N.A,NRHS,NDIMP1.X,NDIMP11.B.IR.IC) 
PARAMETER  (R«8) 


NDIMPl 
C  NDIMPl : 
C  MIAA: 

C  AA: 

C  NRHS 
C  NDIMPl: 
C  DU: 

C  NDIMPl: 
C  RHS: 

C  IR,  IC: 
C 

C  The  input 


«  MIAA  ■  NDIM  <•>!  ■  8  FOR  cab/vaxag  project 
nunber  of  equations,  N 

first  diaension  of  A  froa  DIMENSION  stateaent.M 
N  *  N  matrix  of  coefficients,  A 

nuaber  of  right  hand  sides 
first  diaension  of  U  froa  DIMENSION  statement, N 
on  exit  DU  contains  the  solution  vectorCs),  X 
first  diaension  of  F  froa  DIMENSION  statement, N 
right  hand  side  vectorCs),  B 

integer  vectors  of  diaension  at  least  N 

matrix  A  is  overwritten. 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CSGLE  IMPLICIT  REAL  (A-H,0-Z) 

C 

C  prepares  and  calls  SVDCMP  and  SVDBKSB 
C 

C  SVDCMP  and  SVDBKSB  are  not  included  in  this  code  due  to  copyright  lavs 
C 

C  in  the  remaining  subroutines  U  does  not  iaply  U  the  state  vector  solved 
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c  in  FUlC,  this  U  goss  sith  U,  V,  H  in  using  SVDCMP 
c 

DINEISIOI  Ad.l)  .Ud.M)  .Vfd)  .Vd.H)  .Bd)  .xd) 

DO  20  I-l.l 
DO  10  J>1.N 

10  U(I,J)  -  A(I.J) 

20  COITIHUE 

CALL  SVDCMP(U.N.H.I.H.W.V) 

HNAI-O.ODO 
DO  13  J-l.H 

IFdU)  .GT.HNAX)U1(AX-V(J) 

13  COMTIMUE 
c 

C  *****  sxpsrinsnt  «ith  this  value  belos 

c  sang's  code  does  not  zero  the  small  ej’s  as  advised  in  "lumerical 
C  Recipes  in  Fortran" 
c 

write(6, ♦)*¥■’ , (w(j) , j»l,7) 
tmiH-HMAX* 1.00-20 
DO  14  J>ld 

IFd(J).LT.W!Ii)W(J)»0 

14  CONTIMUE 

CALL  SVBKSBCU.V.V.H.H.H.N.B.X) 

RETURN 

END 

C  . . . 

SUBROUTINE  CAB6(NT0T,LAB,U) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CSGLE  IMPLICIT  REAL  (A-H,0-Z) 

C 

C  writes  u  for  interesting  or  last  points  to  a  file 
C 

DIMENSION  U(7) 

open  (units  1 1 ,  f  ile«>  ’  u .  last .  dat  * ) 
do  10  i»l,6 
write(ll,222)  u(i) 

10  continue 

sritedl,*)  ntot 
222  format(d36.30) 
close(ll) 
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RETURN 

END 


C  . 

c 

C  The  following  subroutines  are  not  used  in  this  ezanple 
C 

SUBROUTINE  BOND (NDIH . PAR . ICP . NBC , UO . U1 . FB . IJAC . DBG) 

C  . 

RETURN 

END 

SUBROUTINE  ICND (NDIM , PAR . ICP , NINT , U , UOLD , UDOT , UPOLD , FI , IJAC , DINT) 

C  . 

RETURN 

END 

SUBROUTINE  FOPT (NDIM , U , ICP , PAR , IJAC . FS , DFDU , DFDP ) 

C  . 

RETURN 

END 
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